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ABSTRACT 



A two-level, quasi-geos trophic atmospheric model on a 
mid-latitude 3"Plane with a single wave in the x-direction 
was coupled with a three level primitive equation ocean 
model through the exchange of sensible heat and momentum 
to investigate the effect of the ocean circulation on tran- 
sients in the atmospheric circulation and vice versa. A 
long-term (79 years) simulation of the atmosphere-ocean 
system was carried out. The experiment was conducted in 
three phases that consisted of integrating the combined. 
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INTRODUCTION 



The use of numerical models to simulate physical con- 
cepts that would otherwise be restricted to limited observa 
tional studies has advanced the knowledge of both the 
atmosphere and ocean and their large scale interaction a 
great deal over the last several years. But that knowledge 
is far from complete and it is the intent of this work to 
reveal more information on the long term, large scale inter 
actions and the frequencies of the response between the 
atmosphere and ocean. The primary goal of this study was 
to identify the effects of the ocean circulation on the 
transient (time-dependent) atmospheric circulation, and 
similarly the effect of the atmospheric circulation on the 
ocean. 

The development of a large scale air-sea interaction 
model that would remain computationally stable throughout 
a long-term numerical integration of the dynamical equa- 
tions of motion was required to achieve this goal. It was 
desired that the method by which the integrations proceeded 
be one which allowed the synchronization of the oceanic and 
atmospheric time steps so that a "one-to-one" relationship 
existed between the two media. Many complex numerical 
models of both the ocean and the atmosphere have been 
developed in the past. However, those models all require 
large amounts of computer time because of the very short 
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time steps required to satisfy the linear computational sta- 
bility requirement. The model developed in this work had 
fewer degrees of freedom in the atmosphere and this resulted 
in faster integration of the atmospheric equations. It also 
utilized certain approximations and assumptions that permit- 
ted the use of longer time increments in the integrations. 
The necessity of employing a technique to simulate time 
synchronization was thereby avoided. 

Many numerical general circulation models have been 
developed in the last two decades. Phillips (1956) made a 
long period forecast with a two-level, quasi-geos trophic 
model, starting with an atmosphere at rest. Diabatic ef- 
fects were included in the equations as a linear function 
of latitude, and such quantities as the intensity of the 
heating, the value of the vertical stability and the type of 
frictional dissipation, were empirically specified. The 
flow patterns which developed were quite realistic. With 
the development of this model, Phillips demonstrated the 
feasibility of numerical simulation of the atmospheric 
general circulation with a simple model. 

Mintz (1958) proposed numerical general circulation 
experiments that specified a heating function in which the 
diabatic heating of the air depended on the air-sea tempera- 
ture difference. He showed that this heating function could 
be used to study the circulation of the atmosphere on a 
planet covered by either land or ocean pr a combination 
thereof . 
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In both of the above works, the partial differential 
equations were formulated in finite differences. Lorenz 
(1960) discussed a method of expressing a dependent variable 
as a series of orthogonal functions. An example of such an 
analysis, in one dimension, is the Fourier series, in which 
the orthogonal functions are sines and cosines. The coeffi- 
cients of these orthogonal functions become the new dependent 
variables and a finite number of variables are retained for 
prediction. Bryan (1959) reported on his work with Lorenz 
in which this technique was applied to the variables in a 
two-layer, quasi-geos trophic model. The variables were ex- 
panded into a few leading terms of a spherical harmonic 
series to conduct a simple general circulation experiment. 

Smagorinsky (1963) reported on an extended period numeri- 
cal integration of a two-level, baroclinic primitive equa- 
tion model which had ”land" as its underlying surface. This 
experiment was made for the simulation and study of the 
dynamics of the atmospheric general circulation. The solu- 
tions corresponding to external gravitational oscillations 
were filtered by requiring the vertically integrated diver- 
gence to vanish identically. This was a major departure 
from earlier work in that many hydrodynamic and kinematic 
constraints were removed. 

Mintz (1964) conducted an experiment in climate simula- 
tion by numerically integrating the primitive equations of 
atmospheric motion, with friction and heating terms included, 
over a long period of time. His two-level model domain 
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encompassed the entire earth, with the underlying oceans 
and continents having different thermal properties. Leith 
(1965) described a six-level numerical model with which an 
attempt was made to compute the hydrodynamic and thermo- 
dynamic evolution of a moist atmosphere on the whole globe. 
The model took into account such external influences as 
solar heating, evaporation and surface friction. 

Smagorinsky, Manabe and Holloway (1965) reported the 
results of a nine-level general circulation model of the 
atmosphere. The primitive equations were used in that 
study in which the atmosphere attained a state of quasi- 
equilibrium. The nine levels of the model were distributed 
to resolve surface boundary layer fluxes as well as radia- 
tive transfer by the prescribed distribution of ozone, 
carbon dioxide and water vapor. This model was an outgrowth 
of Phillips’ two-level quas i-geos trophic model and Smagor- 
insky’s two-level primitive equation model. 

Arakawa (1966) derived a finite difference Jacobian 
expression for the advection term in the difference vorti- 
city equation for two-dimensional incompressible flow that 
properly represented the interaction between grid points. 

He showed that the use of this differencing technique main- 
tained the integral constraints on quadratic quantities of 
physical importance and that this was sufficient to prevent 
nonlinear computational instability. Before Arakawa’s work, 
long-term integrations had been possible only with the use 
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of large friction and diffusion terms. Arakawa’s scheme 
permitted long--term integrations without the requirement for 
such artificialities. 

Kasahara and Washington (1967) described a model of the 
general circulation of the earth’s atmosphere which was 
developed at the National Center for Atmospheric Research 
(NCAR) . A distinguishing feature of this model was that 
the vertical coordinate was height rather than pressure, 
though hydrostatic equilibrium was maintained in the system. 
Experiments with a two-level version of this model were con- 
ducted. Long-term (over 100 days) integrations were conduc- 
ted testing various physical parameters and examining their 
affect on the results. 

Numerical models have also been developed for predicting 
the atmospheric conditions in enough detail to provide 
daily forecasts for given areas. Kesel and Winninghoff 
(1972) reported on the U . S, Navy operational five-level 
primitive equation model that has been used successfully 
for operational forecasts since 1970. This hemispheric 
model perhaps typifies the state of sophistication to which 
numerical weather prediction has progressed. 

Early mathematical models of the ocean circulation were 
very simple as were the early atmospheric models. Munk 
(1950) utilized the oceanic mass transport to examine a 
wind-driven ocean which did not include specification of 
the vertical distribution of density or currents. He ad- 
hered to the linear theory in obtaining analytic solutions 
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to the equations for steady horizontal motions. 

(1954) studied frictionless, steady horizontal flow which 
could occur in a homogeneous ocean of constant depth. The 
effects of the thermohaline structure and compressibility 
were neglected everywhere. By considering this idealized 
case he avoided the analytical difficulties presented by 
the nonlinear terms in the hydr odynamical equations. 

Veronis (1963a, 1963b, 1965, 1966) conducted numerous 
numerical experiments using both a grid point method of 
examination and a Fourier method of analysis and made com- 
parisons of the two methods. He generated flow patterns in 
a rectangular, 3“Plane, barotropic ocean by assuming the 
basin could be divided into an interior region in which the 
dynamical balance was quasi^^eos trophic and a boundary region 
where the dynamical balance was governed by inertial forces 
much the same as Munk . He developed both linear and non- 
linear models and experimentally determined optimum values 
of parameters, such as the Rossby number and the bottom 
friction parameter^ needed for the numerical work. 

Bryan (1963) obtained solutions for a time dependent, 
nonlinear, barotropic model of a wind driven ocean by 
numerical integration. He varied parameters such as the 
wind stress, the Rossby number for the interior flow and 
an effective Reynolds number for the boundary current and 
examined the results. He examined boundary currents and 
the separation of the western boundary current from the 
coast for both linear and nonlinear cases. Bryan and Cox 
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(1967) investigated the ocean circulation using a baro- 
clinic model to numerically integrate the primitive equa- 
tions in almost complete form for different values of 
Reynolds numbers. The near linear solutions for low Rey- 
nolds number were most accessible to a numerical net with 
coarse resolution and gave valuable information for under- 
standing the more complicated regimes at higher Reynolds 
number . 

Gates (1968) integrated the primitive equations for a 
rectangular homogeneous ocean with a free surface on a g- 
plane starting from an initial state of rest. He investiga- 
ted the transient Rossby waves and found a well-marked life 
cycle of amplification and decay as they propagated west- 
ward at about 1 m sec 

Bryan (1969) continued his earlier work with Cox by 
adding new features to their five-level model. He treated 
salinity separately from temperature and allowed sea ice 
formation in the new model. This work is discussed in more 
detail later in the section on the combined model. Haney 
(1971) discussed the formulation of the net downward heat 
flux at the ocean surface, and used this formulation in 
experiments with a six-level baroclinic ocean model similar 
to Bryan’s with emphasis on the equatorial regions. This 
model is representative of the present state of knowledge 
in ocean modeling. 

In recent years research on the problem of atmosphere- 
ocean interaction has progressed to the point where very 



13 



complex atmospheric numerical models have been combined with 
similarly complex ocean models and integrated for extended 
periods of time* The major work in this area combined a 
multilevel, primitive equation atmospheric model including 
the earth’s hydrology reported on by Manabe (1969) with a 
multilevel, primitive equation ocean model (Bryan 1969)* 

The goal of their study, quite different from the present 
one, was to determine the effect of the ocean circulation 
on the time average climate of the atmosphere* Their ex- 
periments produced some notable results such as realistic 
heat transports and reasonable distribution of the pole- 
ward transport of energy* In one sense, however, this ex- 
periment was considered incomplete in that the goal of 
calculating a climate equilibrium with the combined system 
was not completely attained* 

That research was carried out in three steps: 
a* A study of the atmospheric model without the effect 
of ocean circulation* The effect of heat transfer by the 
ocean current was neglected* The ocean had zero heat 
capacity and simply acted as an infinite reservoir of 
moisture for the atmosphere* 

b* A study of the oceanic model without any feedback 
from the atmosphere* The distribution of surface tempera- 
ture, wind stress and precipitation over the ocean were 
given and held constant with respect to time* 
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c. A study of the joint ocean-atmosphere model in which 
the two systems were allowed to interact fully with each 
other . 

A comparison of the results of the three experiments 
yielded some insight into the role of the ocean in maintain- 
ing the climate and of the atmosphere in maintaining the 
thermal and dynamical structure of the ocean. 

In the Manabe-Bryan model more computing time was needed 
to obtain a one day integration with the atmospheric model 
than with the ocean model. Therefore, attempts were made 
to economize computer time by integrating the atmospheric 
equations for only one year while integrating the oceanic 
equations for 100 years. This resulted in a distorted 
response of the ocean circulation to fluctuations of atmos- 
pheric circulation, since an atmospheric disturbance that 
should have an effect for possibly one week influenced ocean 
conditions for up to two years. For example, it took one 
week for a cyclone wave to traverse the oceanic region. 
Because of the synchronization, the oceanic part of the 
model felt that cyclone for as long as two oceanic model 
years . 

Wetherald and Manabe (1972) reported on the effect of 
seasonally varying the solar radiation in the combined 
model just discussed. Many of the results obtained in the 
new version of the model were different, the most signifi- 
cant of which was a warming trend which occurred when 
seasonal variations of the solar radiation were introduced. 
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It should be pointed out that an exact statistical equili- 
brium was not obtained in any part of the region of this 
integration . 

After an extensive review of the literature which in- 
cluded all of the above works, it was decided that to reach 
a state of quasi-equilibrium with a coupled and time- 
synchronous atmosphere-ocean model the simplest atmosphere- 
ocean system which still contained the essential features 
of the dynamics of each media was required. The atmospheric 
model selected was a version of Phillips' (1956) two level, 
quas i-geos trophic model modified by implementing Lorenz's 
(1960) idea of replacing the dependent variable with Fourier 
coefficients. This decision was based primarily on the 
speed with which the equations of this model could be 
integrated. The equations governing the atmosphere could 
be integrated for one model year in about 30 minutes on the 
IBM 360/67. A model similar to the one desired had been 
used in an earlier atmosphere-only study at the Naval Post- 
graduate School by Taylor (1970) and was readily available 
for further modification. The details of the atmospheric 
model are described in Section II. 

The ocean model selected was a modified form of Haney^s 
(1971a) six level baroclinic model. The number of levels 
was reduced to three and the horizontal domain reduced to 
an area in the Northern Hemisphere. The depth was reduced 
and it was further simplified by using the 3-plane approxi- 
mation. While this primitive equation model was fairly 
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complex, it was readily available for modification. De- 
tails of the ocean model are described in Section III. 

The ocean and atmosphere models were joined through 
the exchange of momentum and sensible heat. Unlike the 
model of Manabe (1969a) , the atmosphere contained no mois- 
ture so there was no explicit exchange of moisture. In- 
stead, the effects of evaporation and precipitation 
processes were simulated by other means. These simulations 
as well as the details of the coupling of the two models 
are described in Section IV. 

In a manner somewhat similar to the experiment of 
Manabe and Bryan, the research was conducted in three phases. 
In phase I, equations of the combined model were integrated 
until a state of quasi-equilibrium was achieved. This 
required 79 years of simulation using more than 100 hours on 
the IBM 360/67 computer. The solution of a complex non- 
linear system of equations by numerical integration over a 
long time period as was conducted in this simulation caused 
some concern over the accuracy of the numerical results. 

Error analysis studies by Fisher (1965) and others have 
shown that, even for short time integrations, errors in 
amplitude and phase arise in the numerical solution. This 
is not to say that such errors make the results meaningless 
because even if fairly large errors were present at the 
end of a long-term simulation of the atmospheric circulation, 
the gross features of the solution would not be changed 
significantly. However, the results should be examined more 
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from a qualitative than a quantitative point of view. The 
results in this long-term experiment surely contain numeri- 
cal errors but these results are considered to be satisfac- 
tory for this study of the coupled atmosphere-ocean system. 

One hundred day averages of the atmospheric surface 
stress and surface temperature were calculated and used as 
driving mechanisms for the ocean in phase II, an ocean-only 
experiment. The time of this experiment corresponded to the 
last six years of phase I. Similarly averaged sea surface 
temperatures were used as a forcing function during phase 
III, an atmosphere-only experiment. This also was initial- 
ized so that the time corresponded to the last six years of 
phase I. Details of the procedures used in all three 
phases are described in Section V. 

The resulting data were examined in two ways. The first 
was a synoptic-type study in which maps of both ocean and 
atmospheric parameters were studied and their interactions 
noted. The combined model produced both atmospheric and 
oceanic features that compared favorably with observed fea- 
tures. Energy balance diagrams were constructed for the 
atmosphere and these compared well with existing knowledge 
of the real atmosphere. These synoptic results are shown 
in Section VI A. 

The second method of examining results was a spectral 
method. Time series of data recorded over the last six 
years for selected locations in the ocean and atmosphere 
were spectrally analyzed to determine if any short or 
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intermediate term fluctuations such as those studied by 
Bjerknes (1962) and Namias (1969) appear due to the inter- 
actions. The results of this analysis are discussed in 
Section VI B of this work. By comparing the spectral 
characteristics of the coupled atmosphere with those of 
the atmosphere with fixed ocean surface temperatures, the 
role of the ocean circulation in producing certain tran- 
sient atmospheric responses may be studied. 
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II. THE ATMOSPHERIC MODEL 



A. DESCRIPTION 



The atmospheric model was based on the selective appli- 
cation of the quasi-geos trophic theory to the equations of 
motion. Since the large-scale motions in extra- tropical 
latitudes are quasi-geos trophic in nature, the domain of 
this model was limited to the middle latitudes. The model 
was constructed on a with a central latitude of 45 

degrees. Phillips (1963) discussed geostrophic motions 
that were characterized by approximate uniformity of poten- 
tial vorticity. These motions have characteristic wave- 
lengths of about 4000 km. He used an expansion in the 
Rossby number about latitude 45 degrees to show that the 
3-plane approximation is justified for geostrophic motions 
of this type. 

Phillips (1956) developed a two-level quasi-geos trophic 
model that he considered was the simplest model in which 
heating and friction could be incorporated. The basic model 
described in this thesis followed Phillips’ work in the 
manner in which the quasi-geos trophic equations were derived. 

The quasi-geos trophic vorticity equation on a 3”plane is 

+ ikxvt • + e If - fo If 




9t 

r_JL _ 
^3x 




( 2 . 1 ) 
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where is the component of the stress in the j 4 direction, 

Ty is the component of the stress in the jj direction 
and other symbols have their usual meteorological meaning. 

The model is restricted to two levels with a vertical 
structure as shown in Figure 1. The £itream function and 
are representative of the flow patterns in the upper and 
lower layers of the atmosphere respectively. 
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4 1000 



Figure 1. Vertical structure of atmospheric model. 



Boundary conditions at the top and bottom are: 
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where is a constant drag coefficient and \V is the wind 

at anemometer level. The stress at level four was obtained 
by assuming that \V was a fraction, of the vectorially ex- 
trapolated wind, • Phillips (1956) suggested that for 
simplicity in stress calculations, the anemometer wind speed 
can be assumed to be 70 per cent of the geostrophic wind 
speed \V^. The stress was therefore taken to be 

(2.4) 



where jL = 0.7 

A method for estimating used by Manabe 



(1969a) is 



(2.5) 



o 

where k = Von Karman ’ s constant, 
o 

h = height at which the stress is to be evaluated 
and = roughness parameter. 

The quantity was assumed to be a constant equal to 4.5 
-3 

X 10 . This is larger than values normally used for this 

coefficient, but since depends on where in the surface 

layer the stress is evaluated, the appropriateness of the 
value can best be examined from the resulting momentum ex- 
change between the atmosphere and surface. Values of para- 
meters at level four were obtained by linear extrapolation 
in p. The quantity \V^ was thereby obtained from 



W 



4 




( 2 . 6 ) 



Smagorinsky (1963) derived a formula for calculating 
the Internal stress at level two, by assuming that the 
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stress was proportional to the vertical wind shear. A 
slightly modified form of this formula was used in this 
work and can be written 



= (pk)2 (^)2 ('V3 - W^) , 



(2.7) 



where (pk )2 is an exchange coefficient probably dependent 
on the local Richardson number. An estimate of this from 
Rossby and Montgomery (1935) and used by Smagorinsky was 

50 gm cm ^sec Other values of this coefficient are 

“1 -1 -1 -1 
225 gm cm sec (Palmen, 1955) and 500 gm cm sec (Riehl, 

1951). The choice of 50 gm cm ^sec ^ for this model was 
based on many preliminary numerical experiments with the 
model. The quantity 2 the inverse thickness of the 

750-250 mb layer and it is equal to 0.1266 km If typi- 
cal atmospheric wind values are used, the ratio of to {T^ 

is on the order of ten to one. 

The horizontal domain of the model is shown in Figure 2. 
The east-west extent is 4800 km and the north-south distance 
is 6000 km. The width was chosen as 4800 km to coincide 
with the wavelength of the x-wave that will be discussed 
shortly. The southern boundary is at 18N and the northern 
boundary is at 72N. The use of the 3“Pl^ne over such a 
latitudinal extent may affect the results, but only in a 
quantitative manner . 

If the hydrostatic approximation, the equation of state, 

and the geostrophic condition (tJj = z) were applied to 

o 

the first law of thermodynamics, a quas i-geos trophic 
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version 



9t ''9p^ 



where 



and 



6000 




72®N 



18°N 



Figure 2. Horizontal domain of the 
model. Distance in kilometers. 



of that equation would be written 



+ ik X • V(|^) + |- w 

^ o 




( 2 . 8 ) 



— ] = 2.3 m^sec ^cb 
P 

K = R/Cp . (2.9) 

Q = heating rate per unit mass 

temperature from the standard atmosphere. 
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The above constant value of O, the static stability, was 
used throughout the experiment. 

Equation (2.1) was applied at levels one and three and 
Equation (2.8) was applied at level two. The following 
equations were obtained: 



f 0)^ 
1 o 2 



^ -Mk . • V(V ^ 

^ 3 t 3 t 
3 , y X, 

" “ ® 3p ^3x “ 3y ^1 

(V *3) +<k X V*3 • V(V *3) + B„ 3^ + 



3 t 



3 t 



= _ g i_ r__5L ±] 

® 3p ^3x 3y ^3 



( 2 . 10 ) 



( 2 . 11 ) 
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where A p = 500 
Following Lorenz 



mb . 

(I960) , 




ib^, and \]j„ were defined as follows: 
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(2.13) 



Then Equations (2.10) and (2.11) were added and divided by 
two to yield: 



|-^(V^i|^j^j) +lk X -VCV^ip^^) + Ik X Vip^'- V(V^ 1 |^^) 



+ 



3x 




3x 

• 

3y 



( 2 . 14 ) 
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Equation (2.12) was solved for Equation (2.11) was 

subtracted from Equation (2.10) and the result divided by 
two. The quantity a)^ was replaced and the following was 
obtained : 



|^(V^- + Ik X • V(V^- + Ik X • V(V^i|;^) 



+ e 






M 



o 3x 



2 . ■ 
ILJi Q 
2f ^2 
o 
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A p 



9t 



2 [ 



- [ 



9t 

9x 



9x 

9t 






X . 



9y 



-] 



(2-15) 



2f 



where P 



A pa 



=3.69x10^^ cm^ 



To reduce computing requirements, the model was further 
simplified by expressing the dependent variables, and 

as functions limited to one wave in the zonal direction. 
The wave chosen was of the predominate baroclinic scale. 

Its wavelength was 4800 km, the zonal width of the model 
atmosphere. Any atmospheric model that is developed with 
the intention of reasonably portraying air-sea interaction 
must include baroclinic disturbances containing wave num- 
bers of about five to eight. The baroclinic processes that 
are responsible for the resupply of kinetic energy to the 
atmosphere are reasonably represented by the wave number 
chosen. Phillips (1956) and Smagorinsky (1963) have demon- 
strated in experiments with no forcing functions in the 
east-west direction that much of the disturbance energy is 
isolated in a single predominate wave number. 
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By transforming the system of equations to wave--number 
space in the zonal direction, grid points were eliminated 
in the x direction and the requirement for relaxation 
in order to solve for ip was removed. This type of approach 
had further computational advantages in that no linear 
truncation errors arose in x-space. Also, the fact that 
there were no nonlinear truncation errors in x-space re- 
duced the possibility of nonlinear instability. Therefore, 
errors caused by aliasing due to x-dif f erencing were elim- 
inated. The fields were defined as follows: 

~ A(y,t)cos kx + B(y,t)sin kx (2.16) 

and 

= F(y,t) + C(y,t)cos kx + D(y,t)sin kx (2.17) 

where A, B, C, and D are Fourier amplitudes of the disturb- 
ance, and E and F are the zonal means. The quantity k is 
the zonal wave number and a value of (27r/4800)km ^ was 
chosen . 

In the initial development of the atmospheric model the 
heating was simply proportional to the air-sea temperature 
difference. A sea surface temperature, was specified 

from climatological data and held constant in time. Later 
in the work a more complete form of the heating was used 
and this will be discussed in a later section of this thesis. 

The mean temperature of the layer one-three is propor- 
tional to . By using the geostrophic approximation and 
the equation of state, the following equation was obtained; 



2 7 



2f 

T„ ' = — lit 

2 R 



(2.18) 



where T 2 ' is a departure from the mean temperature at level 
two. To obtain physically reasonable values for the atmos- 
pheric temperature at level two a constant was added to T 2 ' ♦ 



Then 



^2 = ^ ^ 250.3 . 



(2.19) 



By using the standard lapse rate, was extrapolated down- 



ward to obtain the temperature at level four, T^, and 



9T^ 

^4 2 az 



( A z) 



( 2 . 20 ) 



where A z is the standard thickness of the 1000-500 mb 
layer and 3T^/az is the standard lapse rate. Therefore, 
only differs from by a constant. 

Since the static stability was held constant, the heat- 
ing at level two was proportional to the air-sea temperature 
difference and was expressed as 



Qo = "7^ Hi (T - T, ) 

^2 Ap I sea 4 



( 2 . 21 ) 



where is a factor derived by Haney (1971) and is a slow- 
ly varying function of latitude. In this work a value of 
-2 -1 

55 cal cm deg was used and held constant in both latitude 
and time . 

When Equations (2.16) and (2.17) were substituted into 
Equations (2.14) and (2.15) and terms of wave number 2k and 
higher were neglected, equations in the Fourier coefficients 
were obtained. The various terms were separated and 
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coefficients of the cosine terms, sine terms and terms in- 
dependent of X were equated. This yielded six predictive 
equations for A, B, C, D, E, and F. The derivation of these 
equations is contained in Appendix A. 

The equations are: 



3 ,3^A , 2. , , ,3E 3^B , 3F 3^D 

*1 - *‘la 7 TT 37 — 



^ B - D 



3y 



3y 
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3y‘ 
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(2.22) 



(2.23) 
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D(k + y ) 
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dy 



A p 



P 

pC 
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|«^|k^2D-B) - |^[|W4l(2 - |£)3},(2.24) 
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|lvJk2(2D-B)- |^[|Wj(2f^ - |^))J,(2.25) 
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3 ,3^E, _ k 3 , 3^B 3^A 3^D 3^C 

^[— ] - 2 ^ ^7~2 " 



3y ^ 9y 9y 
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(2.26) 
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(2.27) 



where 



9T 

<=1 ■ ' 250.3 - ^ 



- ( A Z) 1 j|- 



B. METHOD OF SOLUTION 

1 . Finite Difference Scheme 

Centered time differences were used for all terms 
except those involving friction and heating. Those terms 
were computed at time (t-At), A forward time step was 
utilized to start the marching process, A finite differ- 
ence scheme developed by Matsuno (1966) was used every 24 
time steps. This was a two-step iteration to simulate the 
backward difference method. This scheme was applied in 
this situation to eliminate the computational mode that 
would arise from using the centered time differencing ex- 
clusively and also helps eliminate nonlinear instability. 

A time step of 30 minutes was used. 
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Since a wave number representation was used in the 
X direction, finite differencing was not required in that 
direction. However, centered space differences were taken 
in the y direction. A uniform grid distance of 200 km was 
used for Ay. 

2 . Boundary Conditions 

Free-slip walls existed at the north and south 
boundaries. The boundary conditions at these walls were 

A, B, C, D = 0 at y = 0, y = w (2.28a) 

and 

= 0 at y = 0, y = w. (2.28b) 

Condition (2.28a) was imposed to force the meridional flow 

to equal zero and correspond to the statement that the 

normal geostrophic velocity vanished at the wall. From 

9u 

the zonally averaged equation of motion, ^ ^ = 0 at the 

o t 

north and south walls. The boundary condition (2.28b) en- 
sures that this was true and was also a requirement for 
proper energy conservation. 

3 . Initial Conditions 

Initial conditions of a mean zonal current with a 
superimposed disturbance were specified. The vertical 
shear and hence the thermal wind were near the critical 
value for baroclinic instability. The initial disturbance 
had a wavelength, previously discussed, and a structure 
such that when superimposed on a baroclinic current with no 
horizontal shear it had the maximum gro\vth rate. 
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The initial conditions were: 



A 

B 

E 



200 

f 

o 



s in 






2 -1 
cm sec 



C = D = 0 , 

1200 — cos C-^) cm^ sec 
TT 



and 

F = 600 ^ cos (-^) cm^ sec ^ » (2.29a) 

TT W 

It was further desired that the mean surface wind vanish 

initially, which can be achieved by requiring that 
— ^^4 3 7 ^ 



or 

E = 2F. 

4 . Time Integration Process 

Different methods for solving the time tendencies 
of the Fourier coefficients were utilized. A Gauss-Jordan 
elimination method was used to solve equations (2.22), 
(2.23), (2.24), (2.25), and (2.27). Richtmyer (1957) dis- 

cussed the use of such a procedure for solving a set of 
equations that had all the elements of its corresponding 
matrix vanish except those on the three main diagonals. 
These equations had that property. The method of solution 
for Equation (2.27) was discussed by Phillips (1956) and 
modified to ensure that the equation satisfied the boundary 
conditions. Equation (2.26) was solved by a direct march- 
ing process utilizing the proper boundary conditions. 
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C • PRELIMINARY EXPERIMENTS 
1 • Hadley Regime 

To examine the effect of the disturbances on the 

model atmosphere, an experiment was conducted which should 

have produced a single Hadley type cell in the meridional 

circulation. A linear sea surface temperature was imposed, 

with warm temperatures in the south. 

The amplitudes of the disturbance were set equal to 

zero (A, B, C, and D) and Equations (2.22) through (2.27) 

reduced to: 

.2 



3y 



/ 9 ^ \ _ 8 r S 

2 ^at'^ By ^ A p 
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ri! _ .2 9F ^ -^r-^ 4(ok) (^) ^ 
^ 2 ^ ^ at 3y^ A p V ^2 3y 



3y 



_g_ 



A p P<=D I <17 - ^l7> ) 

2 r /2f 

y KR 



2f A p 
o 



H. [t - ( 
1 sea y 



R 



(F + C 






(2.31) 



For simplicity the following definitions were made: 
^1 = 



^2 = " aV ^( Pk )2 (^)2 » 

^ _ y^Kg 

and ^ — • 

3 2f A p 
o 



(2.32) 



An assumption that a steady state condition existed 

g 

in the model was made, i.e. *tt’( ) = 0, yielding 

o t 






(2.33) 
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(2.34) 



Since was not equal to zero, 

i'%i - ° If - 



and Equation (2.31) can be written 

2f 
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(2.35) 



or 

2 2 f 

S 1_F _ s H 

2-2 ^31 R 

9 y 
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° F - S-H,T + S-H- C = 0 

3 1 sea 3 1 R 1 



(2.36) 



Since F represents the mean temperature throughout 

the layer, the zonal mean thermal wind can be defined as 

9F 



" 8y 



(2.37) 



Equation (2.36) was solved analytically, using the boundary 



condition that = 0 at y = 0, W. 
The solution is 
U 



F = 



whe re 

U = 



and 
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3T 
c 9y 



w. 



f _ w. 1 

''y 2^ " 



1/2 



1/2 

sinh[(S^S^) (y-^) 
cosh[ (S^S^)^^^ |] 



sea 



^4 ^ , 



2f 






(2.38) 



This solution, along with the condition that A thru 
D equal zero and E equals 2F, was used as an initial state 
to expedite reaching a steady state. The initial mean 



thermal wind obtained from (2.38) was 
U 



= s; 



cosh[(S^S^)^/^(y- |)] 
cosh[ (S^S^)^"'^ |] 



U 



(2.39) 



where = 0.00433 degrees per kilometer and the change in 
the thermal wind with time was examined. 
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The complete equations were integrated in time 
starting from this initial state. After 40 days, values of 
were within two percent of the initial value. The mean 
meridional circulation taken from the numerical solution is 
shown in Figure 3a. It consisted of a single cell with 
northward flow at the upper level and southward flow in the 
lower level. The vertical motion field consisted of upward 
motion in the southern part and downward motion in the north. 
The central portion displayed weak and variable vertical 
motion patterns and were interpreted as being meaningless. 

It follows from the equations, that the quantities A 
through D will remain zero. The method by which vertical 
motions and the'mean meridional motions were calculated is 
described in Appendix B. 

A similar experiment was made with the same linear 
sea surface temperature profile as before but with the 
initial conditions given by Equation (2.29a). This initial 
condition primarily differs from the previous one in that 
it contains a finite amplitude disturbance with wave number 
k in the barotropic part of the motion (A ^ 0) . This in- 
tegration required more computer time to reach a quasi- 
steady state. A three cell structure as shown in Figure 3b 
developed after 230 days. The northward and southward 
motions were stronger and the vertical motion patterns 
were more intense and well defined. 
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La t i f ude 



Figure 3. a. Meridional circulation arrived at in Hadley 
regime experiment. b. Meridional circulation of atmos- 
phere with disturbances included. 

A quantitative comparison was made of velocities 
obtained in the two experiments. Table 1 lists the values 
of the upper level meridional velocity in cm sec ^ and the 
mean zonal velocities in m sec ^ at levels one and three. 
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LATITUDE 


MEAN MERIDIONAL VE- 
LOCITY UPPER LEVEL 




ZONAL 


VELOCITY 


LEVEL .1 1 


LEVEL 


3 


HAD 


DISTURB . 


HAD 


BIST , 


HAD 


DIST . 


72.0 


0 . 0 


0.0 


0.0 


0.0 


0.0 


0 . 0 


68.4 


2.9 


3.4 


20.5 


6.6 


8.0 


-0.1 


64 . 8 


3.5 


12.2 


18.9 


6.0 


6 . 8 


-1.9 


61.2 


3.6 


16 .0 


19 . 2 


5.0 


6.8 


-1 . 1 


57.6 


3.5 


6.0 


19.4 


16.3 


6 . 7 


5.1 


54.0 


3.6 


-5 . 3 


19.2 


27.3 


6.8 


11 . 7 


50.4 


3.6 


-13.7 


19.3 


34.9 


6 . 7 


15.6 


46.8 


3.6 


-22.2 


19.3 


36 . 8 


6 . 8 


16.6 


43.2 


3.6 


-28.8 


19.3 


33.3 


6 . 7 


14.6 


39.6 


3.6 


-23.8 


19.2 


24 . 7 


6.8 


9 . 8 


36.0 


3.5 


-8.6 


19.4 


14.3 


6 . 7 


3.7 


32.4 


3.6 


6 .4 


19 . 2 


7.2 


6.8 


-0 . 3 


28.8 


3.5 


10.1 


19.3 


7 . 7 


6 . 8 


-0 . 8 


'25.2 


3.5 


3.1 


18.9 


13.0 


6 . 8 


0.7 


21.6 


2.9 


1.4 


20.6 


11 . 3 


8 . 0 


1.8 


18.0 


0.0 


0.0 


0 . 0 


0.0 


0 . 0 


0 . 0 



Table !• Comparison values of meridional and zonal velo- 
from Hadley experiment and disturbed atmosphere 
experiment* Positive values are northward or eastward 
flow. Upper level is 500-0 mb layer. Units are ra sec 
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2 . Energy Examination 



Equations for the components of potential and 
tic energy and energy transformations in terms of the 
ier coefficients are as follows; (The derivation of 
equations is shown in Appendix C.) 

2 

/ dy = ^ / F dy , 

/ P dy = / (C^ + D^) dy , 

e 2 

r V ^ 1 f r 2 . /SF. 2, , 
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kine- 

Four- 
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(2.40) 

(2.41) 

(2.42) 

(2.43) 

(2.44) 
(2.45 

(2.46) 

(2.47) 
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(2. 49) 



(2.51) 



(2.52) 



The relationship of these quantities can be expressed by 
the following equations; 



37 I - Cp + 



f-/ Pdy =C -C 

dt e p e 



|-r / Kdy =C -C. - D.-D 

dt e k ei es 



and 



(2.53) 

(2.54) 

(2.55) 



f Kdy = -C +C, - D.-D 

dt z z k zi zs 



(2.56) 



If Equations (2.53) to (2.56) are summed, the following 
equation is obtained: 



~ f (P+K) dy = G - D,-D -D.-D 
dt ^ z ei es zx zs 



(2.57) 



and shows that the time rate of change of total energy in 
the system is equal to that generated minus that dissipated. 

The quantity represents the conversion of mean 
potential energy to eddy potential energy that is accom- 
plished by the horizontal transport of sensible heat by 

the eddies. The temperature in the layer one-three is pro- 

9 F 

portional to and k(BC - AD) is proportional to v'T' 
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dT 

dy' 



The quantity k(BC -- AD) represents v’T’ which is the 



transfer of heat across a latitude circle, and since 



9T 

9y 



tends to be negative, k(BC - AD) must be positive to convert 

P to P . 
z e 

The Term represents the conversion of eddy poten- 
tial energy to eddy kinetic energy by the vertical circula- 
tion in the zonal plane. The conversion depends on the 
correlation between O) ’ and a’ (a* is a departure from a 
mean specific volume of air). Since O) ’ a ’ is proportional 
to 0)^ and C, and 03^^ and D must be negatively corre- 

lated to convert P to K . 

e e 

The quantity represents the conversion of eddy 

kinetic energy to zonal kinetic energy and is dependent on 
the correlation between the mean zonal wind and the meri- 



dional convergence of the eddy momentum transport. The 
g 

quantity u -r — (uW’) is proportional to 

2222 
, B ^3 A . ^3 D C, 

> 



8y 
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9y 



9y 



9y 



9y 



3 E 3 F 3 

Since and tend to be negative, y^(u’v’) must be 

positive in order to have conversion of K to K . 

e z 

The term represents conversion of zonal poten- 
tial energy to zonal kinetic energy by the mean meridional 
circulation. It depends on the correlation between O) and 
a, which is proportional to o) . This conversion is small 
in nature ( Oort, 1964) and requires a positive correlation 

to convert K to P . 

z z 
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The quantity is the generation of zonal potential 
energy through diabatic process. In this model the heating 
is a function of the air-sea temperature difference. The 
correlation of T and Q, or F and Q, must be positive to 
generate zonal potential energy. This means that heating 
of relatively warm air and cooling of relatively cool air 
occurred in the model. 

The dissipation terms D D , D . and D provide 

ei es zx zs 

an energy loss mechanism for eddy and zonal kinetic energies. 
The second subscript, s or i, refers to dissipation by sur- 
face or internal friction respectively. 

An interesting way to examine the system energy is 
by considering an energy flow diagram. Figure 4 depicts 
such a diagram. The boxes contain time rate of change of 
energy, and the arrows show the direction of the flow of 
energy if the transformation is positive. 

The model equations were integrated until a statis- 
tically steady state was reached. Energy calculations were 
made and instantaneous values for the quantities just dis- 
cussed are shown in Figure 5 . 

The diabatic heating from the sea surface creates 

-2 

zonal potential energy that amounts to 1.27 watt m . The 

conversion from zonal potential energy into its eddy 

counterpart occurs at the rate of 1.32 watt m . The eddy 

potential energy is converted to eddy kinetic energy through 

-2 

the eddy circulation at a rate of 1.34 watt m . The 
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Figure 4. Energy flow diagram. Flow is in the direction 
of the arrows if the transformation is positive. 

conversion of eddy kinetic to zonal kinetic energy, which 

maintains the zonal currents against frictional dissipation 
-2 

is 0.63 watt m . The dissipation of eddy kinetic energy 
by surface friction is 0.79 watt m . This is substantially 
larger than dissipation by internal friction and the dissi- 
pation of zonal kinetic energy by surface friction. 

It should be noted that the change of potential 
energy, both zonal and eddy, was very small while the change 
in kinetic energy was substantial. The atmospheric model 
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Figure 5. Flow diagram of the model atmosphere energy. 
Energy and energy transformation units are in watt m“2 
(10^ erg cm* ^ sec*^). 

had reached a statistically steady state at the time the 
energy computations were made, but the energy continued to 
fluctuate about a mean value. This undoubtedly influenced 
these calculations. 

An exact balance was not achieved throughout the 
entire system. It is believed that these very small dis- 
crepancies were due to truncation or roundoff errors and 
should not be considered a serious defect in the calcula- 
tions. Figure 6 is included to show that these values are 
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Figure 6. Flow diagram of atmospheric energy as con- 
structed by Oort(1964). Energy units are in 10^ joule 
m*"^; energy transformation units are in watt m“^ . 



not only qualitatively correct but also quantitatively 
acceptable. Figure 6 shows the atmospheric energy flow 
diagram as constructed by Oort (1964). The values en- 
closed in the boxes are no longer time rates of change of 
energy but the amount of energy present. 

It is satisfying to note that the transformations 
in the model atmosphere compare favorably with those of 
Oort. The conversion and agree very well while 
and lie just outside the extreme values as proposed by 
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Oort. The dissipation of eddy kinetic energy is much small- 
er than Oort’s. Since his value itself is considered too 
small, this points out a weakness in the dissipative mechan- 
ism of the model. This may be due to the lack of horizontal 
diffusion terms in the model equations. It was decided that 
diffusion was unnecessary in the model and this decision may 
be reflected in the energy dissipation. Also, energy loss 
through nonlinear interaction between waves has been elimi- 
nated in this model by the single wave formulation. 
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Ill, THE OCEAN MODEL 



A. DESCRIPTION AND EQUATIONS OF THE MODEL 

The ocean model is similar to that studied by Bryan 
(1969) and Haney (1971a), It is a baroclinic model based 
on the primitive equations in which steady state solutions 
are sought by direct numerical integrations of an initial 
value problem. The circulation produced in the ocean-only 
model is primarily driven by surface forces that simulate 
wind stresses, and horizontal pressure gradient forces due 
to internal stratification. The horizontal stress exerted 
on the ocean surface by large scale wind systems causes the 
surface water to move as an Ekraan layer and drift to the 
right of the stress in the Northern Hemisphere, This causes 
the curl of the wind stress to become important as it causes 
convergence or divergence of the Ekman layer transport and 
corresponding vertical velocities below the Ekman layer. 

Heating in low latitudes and cooling in high latitudes 
causes a thermal circulation to develop, A thermal current 
is set up where there are easterlies in the south and west- 
erlies in the north. These currents converging into the 
continental boundaries contribute to the resulting circula- 
tion , 

The rectangular ocean basin is everywhere 1700 meters 
deep and is 4800 km wide and 6000 km in length. To coin- 
cide with the atmospheric domain the basin extends from 18 
degrees to 72 degrees north (Figure 2) , As in the model 
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atmosphere, the g-plane approximation is used. Veronis 
(1963a, b) has shown that the p-plane approximation is satis- 
factory for the intermediate scale motions in middle lati- 
tudes. This formulation is also used to maintain consistency 
with the atmosphere and prepare for future joining of the 
two models. The 3~plane approximation in both models elimi- 
nates the need for mapping and reduces the possibility of 
confusion that might arise if an ocean model in spherical 
coordinates were joined with an atmospheric model using car- 
tesian geometry. 
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Figure 7. Vertical structure of the ocean model. 



The vertical structure of the three-level model, shown 
in Figure 7, is chai\ged considerably from that described by 



Haney. The three levels are placed at 100, 500, and 1300 
meters below an undisturbed surface. Temperature and hori- 
zontal velocities are predicted at these levels while the 
vertical velocity is diagnostically computed at the inter- 
mediate levels. The vertical layering is defined in such a 
way that there is a finer resolution near the surface than 
near the bottom. 

The Boussinesq approximation (variations of density are 
neglected except in terms involving gravitational accelera- 
tion) and the hydrostatic approximation are applied. In 
the equation of state, the dependence of density on salinity 
is neglected. Density is assumed to be a linear function 
of temperature only. 

The governing equations are 



du _ 1 

d t p 9x 



+ vf + — F 



P X 



(3.1) 



dv _ __ 1__ 0 ^ 

d t p 9y 

o 



uf + — — F 

p. y 



(3.2) 
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(3.4) 



(3.3) 



dT ^ 5 

dt p C 



(3.5) 



and 



p = p^[l - a(T - T^)] 



(3.6) 
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where F = eastward frictional force per unit volume, 

F = northward frictional force per unit volume, 
and Q = heating rate per unit volume. 

The quantities F^ , F^ and Q represent the assumed con- 
tributions of sub-grid scale motions to the exchange of 
momentum and heat in the ocean. The internal redistribu- 
tion is accomplished by utilizing a linear lateral diffusion 
of heat and momentum. Similarly, the vertical mixing of 
these two quantities by sub-grid scale eddies is specified 
by a linear vertical diffusion coefficient. The vertical 
mixing of temperature is enhanced by the use of an instan- 
taneous adjustment mechanism to eliminate unstable lapse 
rates . 

For the initial experiments X7ith the ocean model a con- 
stant wind stress was imposed. Figure 8 shows the mean 
zonal stress that was used in this phase. The data was 
taken from Hellerman (1967). 

Following Haney (1971a) the downward heat flux at the 
ocean surface is expressed as 

- '“o 

The quantity contains the net downward flux of solar 

radiation across the ocean surface, minus the upward flux 

of longwave radiation and latent heat from an ocean surface 

at a temperature equal to T , the prescribed atmospheric 

a 

equilibrium temperature. The second term contains the up- 
ward longwave radiation and sensible and latent heat. The 
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coefficient was found by Haney to vary with latitude by 
about 20 percent. It was set equal to a constant in this 
work . 




Figure 8. Mean zonal stress used in the ocean model 
experiments. The data is from Hellerman (1967); units 
in dynes cm~2 ^ 

B. COMPUTATIONAL PROCEDURES 

Since the ocean bottom is flat, w = 0 is the boundary 
condition on the vertical velocity there. Following the 
method used by Bryan and Cox (1967), Haney (1971a) and 
others, the ocean surface was considered a balanced surface 
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at which w 



0. The height of the surface v/as not obtained 



from the continuity equation so external gravity waves were 
eliminated, which permitted a much larger time step to be 
used. However, the technique also removed the vertical 
mean divergence from all scales of motion, including Rossby 
waves, and introduced an artificial constraint on the phase 
speed of these waves (Gates 1968). The necessity for the 
longer time step to accomplish long-term integrations made 
this technique mandatory. This filtering technique does 
not affect the thermally driven part of the motion. 

A new set of equations was formed by subtracting the 
vertically integrated form of Equations (3.1) and (3.2) from 
Equations (3.1) and (3.2) respectively. After eliminating 
the pressure deviation in favor of the temperature by the 
hydrostatic equation and the equation of state, equations 
for the time rate of change of the vertical shear current 
(u’v’) were obtained. These equations have the form 



Iy(u') = F^(u,v,T) 



(3.8) 



and 




(3.9) 



where 



u'(x,y,z,t) = u(x,y,z,t) - u(x,y,T) 



(3.10) 



and 



v'(x,y,z,t) = v(x,y,z,t) - v(x,y,T), 



(3.11) 
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and u and v are the vertical mean motions discussed shortly. 
Equations (3.8) and (3.9) were numerically solved to obtain 
new values of the vertical shear current. 

A vorticity equation, from which the vertical mean cur- 
rent was predicted, was formed by eliminating the pressure 
between Equations (3.1) and (3.2). The result was vertical- 
ly averaged and had the form 

It 1/ 

-H 



where C is the relative vorticity and H is the depth of the 
basin. Vertical integration of the continuity equation 
(3.4) allowed the introduction of a stream function for the 
vertically integrated velocity since the vertical mean 
divergence was zero. The vertically averaged vorticity 



therefore had the form 

o 



1 / 



Cdz = 



(3.13) 



-H 



and Equation (3.12) can be written 



h = F3(u,v) . 



(3.14) 



A relaxation technique was used to solve Equation (3.14) 
for the time rate of change of ij; . The quantity was then 
set to zero on the boundaries. A time integration was then 
performed to obtain the new value of the stream function 
and then the vertically averaged velocity components from 




and 




(3.15) 
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The total current was then obtained by adding the verti- 
cal shear components predicted by Equations (3.8) and (3.9) 
to the vertically averaged components predicted by (3.14), 

To obtain new values of the temperature, Equation (3.5) 
was written as 

Iy = (u,v,T) (3.16) 

and integrated in time to obtain T. The quantities p and w 
were then diagnostically calculated from (3.6) and (3.4) 

N 

respectively . 

The boundary conditions on momentum and heat at the 
walls were those of zero slip and perfect insulation. The 
ocean bottom was also considered to be perfectly insulated. 
As a simplification of the model, bottom friction was neg- 
lected and free slip was allowed. 

C. HEATING AND FRICTION 

It was assumed that all solar radiation penetrating the 
ocean’s surface was absorbed in the top layer which was 
300 meters thick. Since the side walls and bottom of the 
basin were insulated, the only source of heat was at the 
surface and was calculated from the thermal boundary condi- 
tion (Equation 3.7). A parameterization scheme by which 
the internal redistribution of heat by sub-grid scale 
phenomena was represented by a simple linear eddy diffusion 
of heat and an instantaneous convective adjustment mecha- 
nism was used in the model. The right hand side of Equation 
(3.5) can be written 
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(3.17) 



_2_ = 
P C 



A V^T + 8-(k^) + iilic . 
*H 32^ 3z-^ At 



Bryan and Cox (1967) have shown that the lateral mixing 

coefficient, A , for the apparent temperature should be 
H 

less than that for momentum given by . This is due to 
the larger scales of turbulent motion in the ocean having a 
tendency to take place along isentropic surfaces rather than 
strictly horizontal surfaces. The lateral mixing of momen- 
tum and temperature in the model is somewhat dictated by 
the resolution of the numerical grid. The coefficient 
should be sufficiently small so as not to obscure the 
lateral transfer of the quantities. A minimum value of 
both and A^ is needed to satisfy the zero slip boundary 
conditions, and to simulate the horizontal eddy transport 
of heat and momentum by scales of motion too small to be 
resolved by the numerical grid. In this phase of the work, 
coefficients were held constant with 

A, 



H 



ot; 2 -1 

= . X 10 cm sec 



and 



g 2 2^ 

A., = 1.0 X 10 cm sec 
M 

The vertical diffusion coefficient, K, differs in that 

the exact physical processes by which vertical exchange 

takes place within the thermocline is unknown. Following 

2 -1 

Bryan (1969) a value of 1.5 cm sec was used throughout 
this study . 
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The convective adjustment mechanism (5T)^ is difficult 
to state mathematically. This device is a means by which 
the temperature structure in the model was forced to remain 
stable. After the convective and diffusive heating terms 
were added, the vertical temperature profile was examined 
for unstable lapse rates. If the temperature field resul- 
ted in an unstable lapse rate, the temperatures were in- 
stantaneously adjusted to the vertical mean temperature of 
those unstable layers. When this first adjustment was com- 
pleted, the temperature field was again examined to deter- 
mine if any other layers had been made unstable. The 
adjustment process was repeated until all layers were 

gravitationally stable. 

8 T 

The quantity K-r — was evaluated as follows: 



1 



[H + H, (T 
o l a 



T ) ]/p C 
sea o 



1 



^1 " ^2 
-) 

^1 - "2 



and 



^9z‘'k=2i 



^9z^k=3j 



^2 " ^3 
K(-^ -) 

^2 - ^3 



= 0 . 



( 3 . 18 ) 



The last equality imposed the condition that the ocean 
floor was insulated since the vertical temperature gradient 
vanishes . 

The vertical eddy diffusion of heat at level one, for 
example, was then defined to be 
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dz 



K(T - T ) 

[H + H, (T - T )]/p C ^ 

o l a sea ' Cz," z.) 



^^1/2 ■ ^ly) 



(3.19) 

The heating term due to lateral eddy diffusion of heat in 



(3.16) may be written as 



. n2.„ . r9 /9T- , 3 /3T., 

^ ^ 3y^3y^^* 



(3.20) 



9 T 9 T 

The terms and were defined at points one-half grid 
distance to the east and north, respectively, of points at 
which T was stored (see Figure 9). The boundary conditions 
were 



11 ) 1 . . 11 ) 3 

3x'l- i 3x'l- Y 



and 



il) 1 . _ 11) 1 

9x"^i=i + ^ 9x^i = i - 77 

max 2 max 2 



(3.21) 



where i = — is a point one-half grid distance west of the 

western boundary of the grid(located at i=l) and i + y 

max ^ 

is a point one-half grid distance east of the eastern 

9 T 

boundary (located at i = i ). These conditions force -r — 

m ax o X 

to pass through zero at the eastern and western boundaries 

thereby producing thermal insulation. Similar boundary 
9 T 

condition on t — at the northern and southern boundaries 
3y 

ensure that the normal temperature gradient vanished at 

the walls, making them perfect insulators. 

F F 

X V 

The frictional terms and — ^ in Equations (3.1) and 

P P 

(3.2) are composed of two parts. These are the sub-grid 
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scale vertical eddy stress and the sub-grid scale lateral 
diffusion of momentum. These can be expressed as 



and 



F , 3t 

— = — -r + A^, V u 

p p 3 z M 

o ^o 



F 

_JL = 



1_ 

P, 



3t 



3 z 



+ A„ V^v 



(3.22) 



9 T 

As with the term Kr — , the quantities T and T were defined 

3 z ^ X y 

at the half integer vertical levels. The vertical eddy 
stress was prescribed at the surface and defined to be 
zero at the bottom. The zonal stress was defined at each 
level as 



— ), 1 
Po 2 



prescribed surface stress. 
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Po ^=^2 



= K 
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(3.23) 



and the meridional stress was defined as 



T 

1 

Po 2 



= 0 



T 

Po ^=^2 



= K 



(Vi - vp 



(z- 
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T 

^ 1 = K . V 
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(V2 - V 3 ) 
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(3.2A) 






Thus the prescribed stress was purely zonal. 

The X component of the term involving the lateral eddy 



diffusion of 


momentum can be 


expressed 


as 
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^3x '^9x-' 9y 


(— ) ] 
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(3.25) 


The quantities 
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(3.26) 



These conditions impose the zero slip boundary condition 
on the velocity at the north and south walls and the no 
flux condition at the east and west walls. Similar defini- 
tions for the meridional components impose the zero slip 
condition on the east and west walls and the no flux condi- 
tion on the north and south walls. 



D. FINITE DIFFERENCE TECHNIQUES 

A space differencing scheme that utilized a centered 
differencing with a staggered grid was used in the integra- 
tion of the primitive equations. When a centered space 
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scheme is used with an unstaggered grid, a computational 
mode in space arises that is analogous to the computational 
time mode that occurs with a centered time differencing 
(leap frog) scheme. The use of the staggered grid not only 
prevents the appearance of the computational mode in space 
but is responsible for considerable savings in computing 
time b.ecause variables are retained at alternate grid 
points. Arakawa (1966) has shown that maintaining integral 
constraints on certain quadratic quantities of physical im- 
portance is sufficient to guarantee nonlinear computational 
stability. The space differencing used in this model is a 
simplified, energy conserving version of the scheme used in 
the UCLA general circulation model, as documented by Gates 
e^ al (1971) . 

Total energy is conserved in this model. The boundary 
conditions of zero normal flow through the side walls, and 
zero vertical velocity at the bottom, guaranteed that the 
following integral constraints were satisfied by the differ- 
ential equations. 

Kinetic Energy 

^ y' J v^)dv = J' [-pgw + Pq(uF^ + vF^)]dv. (3.27) 

V V 

Potential Energy 

pg(z+H)dV = J Ipgw - g(z+H) |^]dv , (3.28) 

V V 

The potential energy was taken to be zero at the bottom. 
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If (3.27) is added to (3.28), a conservation equation 
for the total energy in the basin is derived and is given 



by 



[j + pg(z+H)dv = 



V 



V 




# 



(3.29) 



Finite difference forms of the prediction Equations 
(3.8), (3.9), (3.12) and (3.16) are shown in Appendix D. 

Haney (1971a) has shown that these difference equations 
conserve total energy. 

Figure 9 shows the staggered grid network. The velo- 
city components, u and v, are located one-half grid incre- 
ment east and one-half grid increment north of the points 
where corresponding scalar variables T, p and ip are re- 
tained. The grid distance is 200 km in both the x and y 
directions . 



j 




1 2 .... IM-1 IM 



Figure 9. Location of variables in the horizontal 
staggered grid. 
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In the time differencing scheme, various terms in the 
equations of motion were treated differently. The pressure 
gradient term was evaluated using a centered difference. 

The Coriolis terra was evaluated by a trapezoidal implicit 
scheme, and the friction term was evaluated using a forward 
time step to ensure linear computational stability. 

The trapezoidal implicit scheme for the Coriolis acce- 
leration reduces to a form 

- u’^"^ = 2At[aif + (l-a)if (3.30) 

It is easy to show that this scheme is linearly computation- 
ally stable for all A t if a ^ 1. In this model, 

a = 0 .55 . 

To supress computational modes in time that can arise 
when a leap-frog scheme is used exclusively in an oscillat- 
ing system, a Matsuno (Euler-Backward) step was used 
periodically. This two-step scheme damps the computational 
mode and permits the long term integration. 

The leap-frog scheme is conditionally stable and the 

maximum time step was determined by the highest frequency 

governed by this scheme. External gravity waves have been 

removed from the model which leaves internal gravity waves 

as the highest frequency waves present. The phase speed of 

these waves is Vg * H ^ where g’ is a modified gravity equal 

to g — - -P - , In this model g’ ^ 3.60 cm sec H’ = 1200 
^o 

meters and Ay = 200 km, and the linear computational sta- 
bility criterion is 
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.375 day s . 



A t - 

/i^ 

This was the time step used in this phase of the model. 

E. INITIAL CONDITIONS 

A state of rest was imposed on the ocean for the start 
of the first integrations. A stratified temperature struc- 
ture was selected so that a minimum amount of computing time 
would be required to bring the ocean model to thermal equi- 
librium. Figure 10 shows the vertical temperature structure. 




4j0 43 5.1 5.7 63 69 72 7.8 85 







L 


























c 



4.0 40 40 4.0 4.0 40 40 40 40 

1 1 1 1 1 1 1 1 1 



72' 



L a fi t u d e 



18 



Figure 10. a) Initial temperature field at level one; 
b) Same at level two; c) Same at level three. 



62 



The northernmost temperature at each level is 4 degrees 
Centigrade and the vertical structure is isothermal in that 
column. The temperature at level one increased from 4C at 
the northern boundary to 26C at the southern wall. At level 
two the temperature increases from 4C to 8,5C and at level 
three it is a constant 4C throughout. The stability varies 
from isothermal in the northernmost column to very stable 
in the southernmost column. No zonal variations were ini- 
tially imposed , 
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IV- THE COMBINED AIR-SEA MODEL 



A combined model was developed by joining the atmosphere 
and ocean models* The exchange of momentum and heat across 
the air-sea interface is the mechanism upon which the cou- 
pling of the two models was based. An exchange of moisture 
was not possible since the atmospheric model was ’’dry”, but 
a crude attempt at parameterizing latent heat transport in 
the atmosphere was incorporated and will be discussed in 
the section on heat exchange. 

The physical aspects of the combined models remained 
the same as they were in the individual models. The spatial 
grid in the ocean, and time increments for each model, re- 
mained the same. The time marching process allowed the 
time integration to remain synchronized and avoided a dis- 
torted response of the ocean circulation to atmospheric 
fluctuations. The atmospheric solution was evaluated at 
grid points in the x-directlon in the atmospheric model. 

This was done to construct a complete stress field from the 
atmospheric Fourier coefficients in order to calculate the 
air-sea momentum exchange. The ocean surface temperature 
was transformed to wave number space to calculate the air- 
sea heat exchange. The overall techniques of integration 
used in each separate model remained unchanged in the com- 
b ined model . 
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A. TIME SYNCHRONIZATION 



One of the problems that arose in early air--sea models, 
such as that developed by Manabe (1969a, b) and Bryan (1969) 
was the imprac ticali ty of using a marching process that 
would keep the ocean and atmosphere synchronized. Primitive 
equation atmospheric models require time steps on the order 
of minutes and cannot be realistically integrated in sequence 
with ocean models that have time steps on the order of hours, 
because the computing time requirements are prohibitive. 
Techniques have been developed to simulate time synchroniza- 
tion but little success has been achieved in long-term one- 
to-one integrations. 

Since the present atmospheric model is quasi-geos trophic , 
a time step of one-half hour was utilized. The ocean time 
step was nine hours as discussed previously. The resulting 
ratio of atmospheric time steps to ocean time steps was 
eighteen to one. The marching process consisted of 18 
thirty-minute time steps in the atmosphere followed by one 
nine hour time step in the ocean. 

Conservation of momentum and heat in the exchange pro- 
cess during the integration was crucial to the development 
of a realistic interaction model. Therefore, the quantities 
such as the total surface stress and heating calculated dur- 
ing the 18 atmospheric time steps had to equal the quantity 
as it was applied during the one ocean time step. The 
method by v/hich this was accomplished will be discussed in 
the next two sections. 
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B. MOMENTUM EXCHANGE 

The surface stress discussed earlier and computed from 
Equation (2.3) was used as an upper boundary condition for 
the ocean. The calculation of this stress was based on the 
assumption that there was no motion in the ocean. That 
assumption is also used in the stress formulation in the 
combined model. Since the velocity of ocean currents is 
much smaller than the velocity of air flow, the effect of 
ocean currents on momentum exchange can be neglected. 

Manabe (1969b) used the same argument in justifying his use 
of a similar formulation of the stress. Simplicity and 
computing requirements also entered into the decision to 
use this formulation here. 

Starting from a prescribed initial state, the atmos- 
pheric equations were integrated over 18 time steps. The 
surface stress was calculated each time and an average was 
computed to be used for the ocean integration. A general 
numerical procedure for carrying out the computation of 
the time mean surface stress is given by 

^ NM 

STXX.. = , (t )ijn] (4.1) 

ij NM n = l X 

and NM 



where NM = 

and i and j 
directions . 



(At) /( A t) ^ 

ocean atmosphere 



are horizontal grid indices having their usual 
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Since the surface stress components in the atmospheric 



model, T and T , were not calculated directly, terras had 

^ X y 

to be collected and values for the stresses recovered at 
each grid point. The stresses were computed in three equa- 
tions, (2,22), (2,23) and (2,26), 

The following definitions were made: 

T ),, T,,) , = Contributions of the A disturbance to the 

X A Y A 

zonal and meridional components of stress 

computed in the A equation. 

) T) > ^ Contributions of the B disturbance to the 

X B Y B 

zonal and meridional components of stress 
computed in the B equation, 

'T ) r. 9 T^,)^ = Contributions of the mean field to the 

X E YE 

zonal and meridional components of stress 
computed in the E equation. 

To reconstruct zonal and meridional stress fields from 
these quantities, the terms must be multiplied by trigono- 
metric functions as follows: 



(t ).. = T )_ + T ). cos kx . + T )_ sin kx . 

xij xE xA 1 xB 1 



(4,3) 



and 



(t ),. = T )^ + T ). cos kx. + T )„ sin kx . • (4,4) 

y"ij y E y^A r y^ B i 

The locations at which the atmospheric surface stresses 
were calculated coincide with the points in the ocean at 
which the temperature and stream function were calculated. 
Since u and v were calculated at the staggered or half-grid 
points in the ocean system (see Figure 11), it was necessary 



to average the stresses calculated in (4.3) and (4.4) before 
applying the stress to the ocean. The result is 




X i+l,j 



+ (t_) 



X i,j+l 



+ (T^) 






] 



and 



(4.5) 






(4-6) 



STXX and STXY were then calculated from (t )^ and (x )^ 

X ij y ij 



for application to the ocean. This procedure of averaging 
over the number of atmospheric integrations, and applying 
the stress at the proper position was continued throughout 
the duration of the long term integrations. 

C. HEAT EXCHANGE 

To satisfactorily conserve heat in the air-sea exchange 
process, a time averaging procedure similar to that used in 
the momentum exchange was utilized. The equations for the 
heat flux in both atmosphere and ocean must guarantee that 
the total heat lost by the ocean equals that gained by the 
atmosphere . 

The atmospheric heating was previously discussed and 
given by Equation (2.21). In the coupled model, that for- 
mulation was expanded to include the net radiative heating 
of the atmosphere and a crude specification of the heat 
released through precipitation processes. The heat flux 
into the atmosphere was given by 




(4.7) 
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Figure 11. Grid alignment of atmosphere and ocean to show 
non-coincidence of stress points in atmosphere and velo- 
city storage points in the ocean. 
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where R 5 net flux of radiation into the atmosphere, 
and L • P = net flux of heat into the atmosphere due to 
precipitation. The specification of R and L • P are des- 
cribed below. 

The ocean heat flux was also modified from the form 
used in the ocean only model. The necessity for modifica- 
tion arose from the fact that in the ocean-only model (like 
the atmosphere-only model) there was no requirement for an 
exact accounting of heat exchange. The ocean could give 
up any amount of heat into the non-existent atmosphere with 
no adverse effects. Similarly, when heat was required, 
there was an endless supply available from the same non- 
existent atmosphere. In the combined model, the loss of 
heat from the ocean by evaporation was a specified function 
of y only. The quantity in Equation (3.7) was divided 
into a term representing the net downward flux of short 
and long wave radiation, S^, minus the upward flux of 
latent heat, L*E. The equation for downward oceanic heat 
flux, Q , was therefore written as 

Q = S. - L-E - H, (T - T, ) (4.8) 

o N 1 sea 4 

where S._ E net downward radiative heat flux 
N 

and L*E E net downward latent heat flux. 

Figure 12 gives a schematic representation of the heat flux 
into the combined model and the calculated exchange 

(T^^^~T^ ) . The arrows originating or terminating out- 
side the model atmosphere or ocean indicate externally 
specified heat fluxes. 
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Figure 12. Schematic representation of simplified heat 
flux into combined model. See text for definition of 
symbols . 

Annual mean values were used in the radiation terms in 
both Equations (4.7) and (4.8). Wetherald and Manabe 
(1972) recently reported on the effects of seasonal varia- 
tion of the solar radiation on the joint ocean-atmosphere 
model developed by Manabe and Bryan. Their model atmos- 
phere showed a marked warming in high latitude. Although 
a similar result might be expected in this model, varying 
the solar radiation should be the subject of a separate 
investigation and will not be included in this work. 
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From Equations (4.7) and (4.8), the total heat flux into 

the ocean-atmosphere system, Q + Q , is given by 

a o 

Q + Q = S,, + R - L-E + L*P. (4.9) 

When Equation (4.9) was integrated over the horizontal do- 
main of the model, the total prescribed heat flux into the 
combined system was obtained. Time integration of the 
model cannot approach a steady state unless this total pre- 
scribed heat flux into the system is zero. Also, the atmos- 
pheric circulation is insensitive to the value of the total 
heat flux because only an input of available potential 
energy impulses the circulation. Since each of the heating 
components on the right hand side of (4.9) represent essen- 
tially different physical effects, each term was prescribed 
to give a zero net heat flux when integrated over the hori- 
zontal domain. The meridional mean was removed from R and 
and the quantities L*P and L*E were made sinusoidal 
about zero in the y direction. Observed data shows that 
the ocean receives radiation and the atmosphere loses it 
at all latitudes. Removing the mean and using the new 
quantities R’ and ’ simulates the role that the media 
play as both a source and sink of radiative heating. Re- 
moving the mean from both L*P and L*E was necessary to 
insure that the net prescribed heating in both the atmos- 
phere (given by (4.7)) and ocean (given by (4.8)) was 
zero. Thus, there is no net prescribed heat exchange be- 
tween the two media - all of the heat exchange is calcu- 
lated internally from the last term in (4.7) and (4.8) and 
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depends on the air^sea temperature difference. Denoting 
new quantities that have a zero meridional average by 
primes, Equations (4.7) and (4.8) become 

= R' + L*P'+ 11^ - T^) (4.10) 

and 

Q = S ’ - L-E’ - H. (T - T,). (4.11) 

^o n 1 sea 4 ^ 

These terms will be discussed in the next sections. 

1 . Atmospheric Radiation, and Solar Radiation, 

The observed annual mean net atmospheric radiative 
cooling plotted as a function of latitude is shown by the 
dashed line in Figure 13. The data is from London (1957). 
It shows a meridional average loss of heat due to radiative 
effects of 184 langleys day The do t-and-dashed line in 

this figure portrays the deviations from the mean radiation 
As can be readily seen, the deviations from the mean are 
very small and are in all cases less than 10 langleys day ^ 
These deviations were considered too small to be signifi- 
cant and were neglected. Therefore, R* E 0 for future 
calculat ions . 

Figure 13 also shows the observed annual mean solar 
radiation minus the mean upward longwave radiation at the 
earth’s surface, . This quantity with its mean removed, 
S^’, is shown by the dashed line in Figure 14. The values 
of this quantity used in the model are shown by the solid 
line in this figure. A profile that was both linear aiid 
symmetric was chosen and held constant in time. 
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Figure 14, Deviation from mean net radiation, (dashed 

line) as observed by London (1957); model annual mean net 
radiation, (solid line). 
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2 , 



Latent Heat Terms, and 



The parameterization of the latent heat processes 
could be accomplished in one of three ways. Since the at- 
mospheric model contained no moisture, a direct exchange 
of this quantity and the associated releases of heat could 
not be considered. One possible formulation was based on 
the assumption that the air directly above the sea surface 
was saturated so no evaporation existed (L*E’= 0), This 
formulation coupled with the assumption that no precipita- 
tion existed thereby eliminated the terms L*E’ and L*P’ 
from Equations (4.10) and (4.11). An experimental short 
term integration of the model equations was completed using 
this formulation. Results were unsatisfactory in that it 
became evident in this short test that an active air-sea 
model could not be developed with this scheme. The merid- 
ional gradient of the sea surface temperature became weak 
and the atmosphere consequently became bar oclinically 
stable. After a short time this formulation was discarded. 

A formulation of L^'E’ and that was similar to 

the observed profiles (Figure 15) was used in a second 

short term experiment. The result was similar to the case 

just described, only the meridional temperature gradient 

was weakened at a more rapid rate. The solar radiation 

absorbed at the ocean surface, ’ was insufficient to main- 

N 

tain an adequate meridional sea surface temperature gradi- 
ent and the addition of latent heating in the north and 
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cooling in the south further reduced the atmosphere’s 
meridional temperature gradient. The results of this test 
were also unsatisfactory. One suggestion for the reason 
that the ocean did not maintain a strong temperature gradi- 
ent even though the applied solar radiation was close to 
that observed in the real atmosphere is that the model can- 
not simulate the sharp separation of the western boundary 
current from the coast and the resultant tight gradient of 
temperature in the mid-latitude regions. This was probably 
due to the degree of linearity of the solution in this area 
(see Section VI A 1) and resulted in an oceanic and atmos- 
pheric meridional temperature gradient that was too weak 
for baroclinic instability. The system as a whole needs 
the latent heat released in the tropics of the real atmos- 
phere . 

The results of the above two experiments indicated 
that the ocean-atmosphere model required a somewhat differ- 
ent partitioning of the heating than is typically observed 
in the real ocean-atmosphere system. The heating formula- 
tion adopted for the model was equivalent to assuming that 
evaporation occurred in the ocean part of the model and 
that the moisture immediately precipitated out at the same 
latitude in the atmosphere (L*E’ = L*P’). This formulation 

is somewhat unrealistic since, over the latitudinal domain 
of this model, the observed evaporation exceeds precipita- 
tion in lower latitudes while the reverse is true in the 
middle latitudes. Figure 15 shows curves of the latent 
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heat flux due to evaporation (solid line), and precipita- 
tion (dashed line), from Budyko (1962). The meridional 
average values have been removed from these values. The 
sinusoidal form of both L*E’ and L*P* that was used to speci- 
fy these fluxes in the model is shown by the do t-and-dashed 
line in Figure 15. While this formulation is quantitatively 
in agreement with the evaporation curve, it is a poor ap- 
proximation to the precipitation related heat release. How- 
ever, the greatest precipitation is in the equatorial region 
lying south of the wall at 18N. If the hemispheric preci- 
pitation were mapped into the domain of the model, the pre- 
scribed latent heat flux due to precipitation would be more 
in agreement with observed values. This is considered to 
be a weakness in the heating formulation and is an area 
that could be improved on in future work with this model. 
Nevertheless, it was believed that once the heating gradi- 
ent in the atmosphere was sufficient to maintain it in a 
baroclinically unstable state, the resulting air-sea inter- 
actions would be realistic. That is, the interactions 
depend on the characteristic atmospheric motions, independ- 
ent of the external forcing that maintain those motions. 

Figure 16 shows the net prescribed atmospheric heat 
flux, R’ + L*P* (dashed line) and the net prescribed ocean 
heat flux, - L*E (solid line). Prior to its final 

acceptance this form was tested in a short term experiment 
and found to satisfactorily contribute to the conditions 
desired for an active air-sea model. This fact coupled 
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Figure 15. Latent heat flux due to evaporation (solid 
line); latent heat flux due . to precipitation (dashed line) 
as observed by Budyko (1962), >Iodel values of these 
quantities (do t-and-dashed line). All have the meridional 
averages removed. 
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Figure 16. Net prescribed atmosphere (dashed line) and 
ocean (solid line) heat fluxes. 
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with the requirement for simplicity led to its use in the 
remainder of this work. 

3. The Sensible Heat Transfer, Ht (T - T,) 

The quantity " *^ 4 ) will be referred to as 

the sensible heat. In the derivation of by Haney (1971), 
it was shown that this quantity actually contained the net 
upward flux of longwave radiation and sensible and latent 
heat which occurs due to the air-*sea temperature difference. 
Since the parameter H^ varied little, a constant value was . 
assumed. As in the ocean-only case a value of 55 langleys 
day ^ was chosen. 

As shown above, all of the heat exchange between 
the atmosphere and the ocean must be accomplished through 
the sensible heat term. It is in this terra that a strict 
accounting of heat exchange had to be done. 

The sensible heat flux applied during an ocean time 
step was the average of the heat flux calculated during the 
18 previous atmospheric time steps. The equation utilized 
was 



NM 



HEAT. . = 
iJ 



H, (T - T. ) . . ] , 

NM t 1 sea 4 ijn ’ 
n = l ^ 



(4.12) 

where NM = (At) ocean/ (A t) atmos = 18. Because of the 
x-wave formulation used in the atmospheric model the quan- 
tity (T - T,) had to be calculated in several steps. 

By using Equations (2.17), (2.19) and (2.20) an equation 



for T^ could be written; 



(T.).. = (T,F) +(T,C) cos kx. +(T,D) sin kx . , (4.13) 

4ij ^4j X 
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(4.14) 



To be consistent with the wave number representa- 
tion in the atmosphere and to facilitate coupling the 
models through a single wave number, a Fourier analysis was 
performed on the sea surface temperature. At each latitude, 
the sea surface temperature was analyzed in x-space to 
yield a mean zonal sea surface temperature and the Fourier 
coefficients of its first mode. This can be symbolized as 

(T ).. = (T ). +(T C).cos kx.+(T D).sin kx..(4.15) 
sea ij sea j sea j i sea j i ^ ^ 

The various modes of the sensible heat flux in the atmos- 
phere were calculated from 

H, (T - T. ) = H, [ (T - T.F) . 

Isea 4 1 sea 4 j 

+ (T C - T, C) . cos kx . 
sea 4 J 1 



+ (T D - T.D). sin kx.] 
sea 4 J 1 



(4.16) 



and were applied in Equations (2.24), (2.25) and (2.27) 

during each time step. The right hand side of Equation 
(4.16) was averaged over the 18 atmospheric time steps for 
incorporation in Equation (4.12). This in turn was applied 
in the oceanic heating formulation. 

The model equations were integrated for a short 
time period (300 days) and zonally averaged values of the 
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sensible heat were calculated. The results are shown in 
Figure 17 by the solid line. The dashed line shows the ob- 
served zonal mean sensible heat flux as presented by Budyko 
(1962). These were added to the imposed values of L*P’ to 
obtain the zonally averaged net atmospheric heating. Values 
of the observed sensible heat flux estimated by Budyko were 
added to the observed values of the radiative and latent 
heat contributions, and after their meridional averages 
were removed, these new quantities were plotted in Figure 
18. The solid line shows the observed atmospheric heating 
over the ocean and the dashed line shows observed atmospher- 
ic heating over ocean and land combined. The dot-and- 

dashed line depicts corresponding values of L*P’ + H- (T 

1 sea 

- T^) generated in the model. In the northernmost segment 
of the domain, the model closely resembles the observed 
values over the ocean. From a region near 45N to the 
southern boundary, the model is in closer agreement with 
the data from a land and ocean environment. 

There are several possible explanations of the 
**dip” in the model atmosphere heat flux near 50N. The 
formulation of the net solar radiation as shown in Figure 
13 favors a sea surface temperature structure that would 
be colder than the atmosphere in the higher latitudes and 
warmer than the atmosphere in lower latitudes. This would 
lead to negative contributions from the sensible heat term 
in higher latitudes. The curve dips sharply in the area 
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Figure 17. Heat flux due to sensible heat exchange. 
Positive is flux from the ocean; Model (solid line); 
Observed by Budyko (dashed line). 
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Figure 18. Observed zonally averaged atmospheric heat flux 
over ocean (solid line); observed zonally averaged atmos- 
pheric heat flux over land and ocean, combined (dashed line) 
model zonally averaged atmospheric heat flux (dot-and-dashed 
line); (observed data from Budyko , 1962). 
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where the solar radiation and sensible heat passes from 
positive to negative values. 

The formulation of the latent heat released by pre- 
cipitation in the atmospheric model contributes to heating 
in the south and cooling in the north. The sharp dip in 
the total atmospheric heating occurs at the latitude where 
the prescribed atmospheric heating becomes negative. 

The sensible heat term may contribute to this 
abrupt decrease in another way. Large positive values of 
this quantity are seen in some regions and large negative 
values in others. The large negative values are physically 
unreal because the real atmosphere does not give large 
amounts of heat to the ocean through this term over large 
areas. These values are considered departures from a mean 
and can be interpreted as small heating where the values 
are negative and larger heating where they are positive. 

The prescribed heating of the combined model has already 
been shown to meridionally integrate to zero. There can 
be no net heat flux by the sensible heat term if a steady 
state condition is to be reached. If a different formula- 
tion for the prescribed quantities was chosen so that the 
net specified heat flux into the system v;as not zero, then 
the net heat flux by H, (T ^ T , ) would have to be adjus- 

ted. In the present formulation, however, there is no need 
to change the mean temperature of the atmosphere and the 
interpretation just discussed is justified. 
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Figure 19 shows the result of subtracting the sen^- 
sible heat term from prescribed radiation and latent heat 
contributions in the ocean model. The curve has the 
general appearance of the curve in Figure 14, But a 

”hump** appears in the same area as the "dip" in the atmos- 
phere curve and should be expected for the same reasons. 

The oceanic heating formulation that was arrived 
at empirically and through experimentation was utilized in 
the remainder of the numerical study. 
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Figure 19. Average of heating in the ocean model. 
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V. 



PROCEDURES 



A. PHASE I 

1 , Variations of the Combined Model 

It became evident during the short term experiments 
discussed in the last section that the computing time re- 
quired to achieve a state of quasi-equilibrium was excessive. 
Bryan (1969) has shown through scale analysis that the ver- 
tical advective time scale is on the order of decades while 
the vertical diffusive time scale is on the order of cen- 
turies. Both of these time scales must be relevant in the 
integration of the coupled ocean and atmosphere over long 
periods of time. The relatively shallow ocean depth and 
stratified temperature structure specified as the initial 
conditions enhanced the approach to an equilibrium state. 

The combined model with a 200 km grid distance and 
time increments as previously described required approxi- 
mately one hour of computing time for every 100 days of 
simulation on the IBM 360/67 computer. This would require 
365 hours of operating time for 100 years of simulation in 
this phase of the experiment, and this was prohibitive. 

Several ways to reduce computing time requirements 
were considered. These considerations involved the ocean 
model since the atmospheric model was formulated simply 
and could be integrated rapidly. In a combined model such 
as this, almost any change in either time increment or grid 
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distance in one part of the model would require an adjust- 
ment of that quantity in the other portion of the model. 

A time reduction could be realized by linearizing 
the horizontal equations of motion, as was done by Haney in 
the first phase of his ocean experiments. This change was 
incorporated into the model, along with the assumption that 
the vertical mean motion was steady so there was no require 
ment to relax for the time rate of change of the stream 

function. It was this feature that produced the time reduc 
tion. This simplified form of the combined model was 
tested in a several hundred day simulation and the model 
advanced 140 days after one hour of computing time. While 
this was a fair improvement in computing requirements it 
was not enough to warrant changing the model in this way. 

Changing the distance between grid points was con- 
sidered. The model previously described had a 200 km grid 
distance in both x and y requiring horizontal fields of 775 
points (25 by 31). If the grid distance was changed to 
400 km in the y direction and remained 200 km in the x 
direction, the computing time would be reduced. The reduc- 
tion was due to the fewer number of calculations to be made 
(now 400 grid points) and the doubling of the time incre- 
ment in the atmospheric portion of the model. For predic- 
tion purposes the atmosphere had grid points only in the y 
direction and if this increment was doubled, the time step 
could also be doubled and the solution still remain stable 
in a linear sense. The time step in the ocean, however. 
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must remain the same since the grid distance in x was not 
changed. This version of the model was tested and produced 
a 320 day forecast in one hour. This was a substantial 
time savings but it was believed that an even more acceler- 
ated pace could be obtained. 

The final modification to the original combined 
model was made by changing the grid distance in both the x 
and y directions to 400 km. When the time steps were 
doubled to one hour in the atmosphere and 18 hours in the 
ocean, computing time was reduced by a factor of about 
seven. A test showed that this version of the model pro- 
duced a 666 day simulation in one hour of computing time. 
This was a satisfactory improvement and this form of the 
model was utilized in the early stages of the experiment. 

2 . The Coarse Grid (400 km) Model 

Munk (1950) has shown that the development of a 
western (frictional) boundary current of several hundred 
kilometers is to be expected. It was assumed that the 
western boundary current was the width of the new grid 
distance, 400 km. It has been shown by Munk that the width 
of the boundary current (L^) is proportional to the cube 
root of the ratio of the lateral momentum diffusion coeffi- 
cient to the gradient of Coriolis parameter. This can be 
expressed as 

\ cC • (5.1) 
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L^,/Ax, is proportional to I 3 ( A x) ^ . In order 



W 



The ratio of the western boundary v;idth to the grid size, 

to 

keep this ratio the same in both the 400 km and 200 km grid 
models, and thereby have similar horizontal resolution in 
both models, the following must hold, 



M 



= 8 A 



(400) 



M 



( 200 ) 



(5.2) 



8 2 — 1 

Since a value of 1 x 10 cm sec was used for A , with the 

M 

8 

200 km grid, = 8 x 10 was used as an initial value in 

M 

this model. Similarly A^^ was increased so that a value of 
8 2 -1 

A = 1 X 10 cm sec was used in the coarse grid model. 

H 

The initial conditions discussed previously were 
modified to conform to the new 13 by 16 point grid and the 
long term integration was started. Because of the large 
horizontal eddy diffusion of heat, the meridional sea sur- 
face temperature gradient dropped from its initial value of 
22 degrees over the 6000 km length to 18 degrees in the 
first 400 days. To have an active air-sea interaction 
model the meridional temperature gradient in the ocean must 
be large enough to force strong meridional temperature 
gradients in the atmosphere so that a thermal wind larger 
than that required for baroclinic instability can exist. 

Such a large temperature gradient could be maintained in 
the ocean only by reducing the diffusion coefficient. How- 
ever, that coefficient cannot be arbitrarily reduced because 
in the reduction, the western boundary width as given by 
(5.1) may become smaller than the grid size. 
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A stepwise reduction in the lateral diffusion co- 
efficient was made* The solution was examined for computa- 
tional stability after each step. At the end of the first 
400 days, the lateral diffusion coefficients were halved, 
the decrease in meridional temperature gradient ceased and 
the solution remained stable* At the end of 1200 days the 
coefficients were again reduced and at the end of 7200 days 

g 

they were reduced a final time to values of A^, = 1 x 10 

M 

g 

and A = .25 x 10 with no effect on the computational 
H 

Stability. It was evident, however, that with the 400 km 
grid, the western boundary current was only marginally re- 
solved . 

The model equations were integrated over the coarse 
grid for a period of 20,200 days (^-55 years). The volume 
average of temperature in the ocean basin was one parameter 
used to measure the approach to equilibrium. This water 
temperature had slowly cooled from an initial value of 
6.3076 degrees to 5.2017 degrees at day 20,200. In the 
early stage of adjustment the cooling was at a very rapid 
rate, 5.2 degrees per century, but by the end of this part 
of the experiment, the net cooling rate had reached zero. 

Another parameter monitored was the meridional 
gradient of sea surface temperature. This quantity had 
stabilized at 21.6 degrees per 6000 km after day 15,000 
and remained so until the end of this time period. During 
this time the atmosphere remained very active. A thermal 
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wind, (U^ -- U^/2) , with maximum values greater than 11 m 
sec which is well above the critical value for baroclinic 
instability, was present in the mid-portion of the domain. 

A strong upper level jet and weaker lower level jet were 
present at the same latitude. The state of the atmosphere 
and ocean was very encouraging and the data for day 20,200 
was used to initialize the fine grid (200 km) model to con- 
tinue the experiment. More discussion of the state of the 
atmosphere-ocean system at that time is included in the 
section on discussion of results. 

3 . Data Retention 

Prior to continuing with the description of the pro- 
cedure used with the fine grid model, a brief discussion 
of data retention must be included. In a project requiring 
large amounts of computing time, it is mandatory that a 
great deal of consideration be given to which data should 
be saved. One must decide not only what parameters to 
retain, but how often and in what form. It is virtually 
impossible to repeat an experiment to extract parameters 
that were not saved. 

The data obviously must be stored for analysis at 
the end of the experiment. Data must also be retained for 
use in restarting the model. It should not be expected 
that a long-term integration of this nature could be com- 
pleted in one computing time period. The equations must 
be integrated and data stored periodically so that after a 
reasonable operating time, such as one hour, the model can 
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be stopped and restarted at a later time using the last 
data in storage as initial values. 

Throughout the experiment with the coarse grid 
model, data were stored at 100 day intervals. The ocean 
variables saved were the stream function from which the 
vertical mean current was obtained, the zonal vertical 
shear current, u’, the meridional vertical shear current, 
v’, and the temperature, T, at all three levels. The at- 
mospheric variables saved were the Fourier coefficients, 

A through F. The retention of such a large quantity of 
data posed a problem, but it was remedied through the judi- 
cious use of magnetic tapes. As one section of the experi- 
ment was completed, the data that would not be needed until 
the final analysis was moved from disk pack temporary 
storage to a set of magnetic tapes where it was permanently 
stored. The necessity for this became even more evident in 
the fine grid model integrations when the fields were in- 
creased from 208 data points per level to 775 points per 
level. Saving this much data would have been impossible 
without a large storage capability. 

The last problem concerning data storage arose dur- 
ing the final part of the long term integration. After a 
state of quasi-equilibrium was approached, five points were 
chosen for retention of daily data for the last six years 
of the experiment and were carefully examined by spectral 
techniques. Data from an ocean-only and an atmosphere-only 
experiment were also saved over the last six years. 
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4 . Fine Grid (2QQ km) Model 

The data for day 20,200 from the coarse grid model 
was linearly interpolated in both x and y to obtain initial 
fields for the fine grid model. The lateral diffusion co- 
efficients for both momentum and heat were held at the 
final values as used in the coarse grid model. While these 
coefficients might have been reduced over a time period as 
they were in the coarse grid, this was not done in order to 
avoid the readjustment that would be required. A further 
change in the diffusion coefficient would increase the ad- 
justment period. It has been shown by Manabe (1969a), 

Gates (1968) and others that the values chosen previously 
are appropriate for the grid size used. 

The equations were integrated over the fine grid 
model for another 6500 days with data retained every 100 
days. At that time, 73 years, the ocean-atmosphere model 
was approaching an equilibrium state. The volume average 
of the temperature initially started to slowly decrease 
but after a few hundred days it commenced an increasing 
trend that it maintained throughout the remainder of the 
experiment. The warming trend began at the rate of 0.8 
degrees per century but this decreased to about 0.5 degrees 
per century at 73 years. Compared to the warming rate of 
2 degrees per century arrived at by the ocean in Manabe 
and Bryants combined model, this was considered an accept- 
able time to start the last six years of the combined 
model experiment. 
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Beginning at the 73rd year, a final 2200 day (6 
yeara) simulation was produced with the fine grid model. 

It was during this time that not only were all the data 
recorded every 100 days but for five selected points, the 
data were recorded every 1.125 days (every three oceanic 
time steps). These points were selected as being represen- 
tative of features that existed in various parts of the 
ocean. Figure 20 shows the location of the five points. 




Figure 20. Location of points where selected data were 
recorded every 1.125 days. 
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lettered A through E. Point A is the northernmost point in 
the eastern sector where sea surface temperatures were 
warmer than the initial values and where strong sinking 
motion had occurred. Point B is in the center of the domain 
where the cyclonic scale atmospheric disturbances had their 
greatest effect. Point C is in the western boundary current 
and, if the model domain is likened to the North Atlantic 
Ocean, this point is close enough to an observation site 
used by Woods Hole Oceanographic Institution so that some 
comparisons can be made. Point D is an area in the south 
central portion of the domain and would be equivalent to 
an area near Bermuda. Point E is in the warm core of the 
ocean close to the southwestern corner of the domain. The 
data at these points were examined spectrally and the dis- 
cussion of this technique is included in a later section. 

At the conclusion of this part of the experiment, 
the combined model had been used for a 79 year integration 
of the equations. This completed phase I of this project. 

B. PHASE II 

This phase of the experiment consisted of integrating 
the equations of the ocean-only model for a six-year period 
using the fine grid. To simulate the coincidence of this 
time period with the final six years of the combined model, 
data from the combined model were used for initialization. 
Initial fields of temperature and velocities were taken 
f rom day 2 6,700. 
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A fixed surface stress and an atmospheric temperature 
had to be imposed. Rather than use climatology, the input 
stress and temperature were obtained from the combined 
model. In preparation for this phase of the experiment, 
the surface stress, and atmospheric surface air temperature 
were averaged over a 100 day period centered at day 26,700. 
These average values were then held fixed during the entire 
time integration of the ocean-only model. 

Data were retained in exactly the same manner as in the 
combined model. Full fields of the ocean variables were 
stored every 100 days and data for the five special points 
were stored every 1,125 days. The variables were examined 
and comparisons were made with the results of Phase I, 

These results are discussed in a later part of this work. 

C, PHASE III 

The equations for the atmosphere-only model were inte- 
grated over the same six year time span with the fine grid. 
Initial values of A through F were taken from day 26,700 
of the combined model data and the 100 day average sea sur- 
face temperatures were imposed as a lower boundary condi- 
tion. Data were saved at the same times as in the first 
two phases , 

Calculations of energy and energy transformations were 
made and are described later. The data from the five points 
were analyzed spectrally and this too is presented in a 
later section. With the completion of this phase of the 
experiment, the long-term integrations were concluded. 
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D, REVIEW OF COMPUTING TIME REQUIREMENTS 

Since the integrations were conducted wi 
numerical models, it is germain at this poin 
computing time requirements for each model, 
presented in tabular form in Table 2. 



th a variety of 
t to review the 
This review is 



MODEL 


DESCRIPTION 


LENGTH OF SIMULATION 
FOR 1 HR COMPUTING TIME 


Comb ined 
(Fine Grid) 


Ax=Ay=200km 

At ^ =30 min 
atm 

At =9 hours 

ocn 


101 days 


Combined 
(Coarse Grid) 


Ax=Ay=400km 

At =1 hour 

atm 

At =18 hours 

ocn 


666 days 


Ocean Only 
(Fine Grid) 


Ax = Ay = 2 OOkm 
At= 9 hours 


308 days 


Atmosphere Only 
(Fine Grid) 


Ay= 200km 
At= 30 min 


697 days 


Combined 
(Linear) 
(Not used) 


Ax = Ay = 2 OOkm 

At ^ = 30 min 
atm 

At =9 hours 

ocn 

linear momentum 
equations 


140 days 


Comb ined 
(Mixed Grid) 
(Not used) 


Ax=200km 

Ay=400km 

At ^ =1 hour 
a tra 

At =9 hours 

ocn 


320 days 



TABLE 2 . Review of computing time 
variations of the numerical models 



requirements 
used in this 



f or 

experiment . 
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It can be readily deduced from these time figures that 
the coupling processes required lengthy calculations. The 
method by which the sea surface temperatures are transferred 
into wave number space required many time consuming opera- 
tions for every oceanic time step. The other summations 
and averaging procedures all contributed to a joining pro- 
cess that caused the combined model to run significantly 
slower than either of its individual components. 
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VI. DISCUSSION 



The data that was collected during the long term inte- 
gration can be placed in two categories. The full fields 
of data and the data collected at the five designated 
points differ widely in frequency of observations and area 
of coverage. Because of this, the methods by which the 
data were to be examined differ considerably. The full 
fields of data that were retained at 100 day intervals were 
examined in a synop t ic- typ e study, while the five-point 
data, recorded at daily intervals, were examined using 
spectral techniques. 

In the synoptic study the analyzed data were compared 
not only to observed climatological data, but also to the 
results obtained from the various phases of this numerical 
experiment. The atmospheric energy balance was examined 
and comparisons were made between phase I and phase III and 
with the work discussed in Chapter II. In the spectral 
study, power spectrum estimates that were obtained from data 
at one particular point were compared with estimates ob- 
tained from data at other points during the same phase of 
the experiment. Similarly, spectral estimates of one 
variable were compared with estimates of the same variable 
at the same point in a different phase of the work. Some 
comparisons were made with spectral estimates calculated 
from observed data collected by the Woods Hole Oceanographic 
Institution , 
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A. SYNOPTIC STUDIES 

Throughout the experiment full fields of the stream func- 

f 

tion, vertical shear current components, and temperatures 
were retained as oceanic data and the Fourier coefficients 
were saved as the atmospheric data. These variables, along 
with fields of oceanic vertical motion and atmospheric zonal 
velocity, vertical motion, mean meridional velocity and 
disturbance amplitude and phase information were printed 
every 50 or 100 days. Many of these variables were examined 
in detail and discussion of the results is contained in the 
subsequent sections. 

1 . The 400 km Grid Combined Model Experiment 

The 400 km grid combined model was merely a tool that 
was used to expeditiously integrate the model equations to 
reach a state of quasi-equilibrium. The ocean-atmosphere 
interaction at the end of this part of phase I was not of 
primary interest. The state of the system at the end of 
day 20,200 (55 years) was to be used only to initialize the 

fine grid model from which more realistic and more meaning- 
ful results may be obtained. However, a brief discussion 
of the conditions to which this form of the combined model 
evolved is necessary to explain why the 200 km grid was 
needed to complete this phase. 

The initial conditions for this part of the experi- 
ment have previously been described. The initial adjustment 
of sea surface temperatures was due to the large hox'izontal 
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eddy diffusion coefficient then in use. The reduction of 
this coefficient and the necessity for maintenance of a 
large meridional sea-surface temperature gradient was al- 
ready discussed. The coefficient of lateral diffusion 
could not be made small enough with the 400 km grid so that 
nonlinear terms would be important, Bryan and Cox (1967) 
have defined an "effective" Reynolds number, Re, as the 
Reynolds number associated with the lateral turbulent mix- 
ing, The nonlinearity of the solution is dependent on the 
size of the Reynolds number defined as follows: 



Re 



_ 




( 6 . 1 ) 



where is the velocity component depending on the density 
gradien t , 

L is a scale length, and 

is the coefficient of lateral momentum diffusion. 

M 

Bryan and Cox conducted a series of experiments 
with values of Re = 0,8, which gave a highly diffusive, 
linear solution and Re = 8,0, which yielded a moderately 
diffuse, slightly nonlinear solution. The maximum value 
of Re used in this experiment was 6.5 and it is believed 
that a value greater than 10 was necessary to achieve a 
truly nonlinear solution. It is easily shown that with the 
400 km grid, continued reduction of produces a situation 
where the linear western boundary current becomes unresolv- 
able [see Equation (4.1)] prior to the occurrence of non- 
linear effects. Therefore, the solution with the 400 km 
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grid is similar to the viscous linear solution described 
by Munk (1950), Figure 21 shows the ocean stream function 
obtained at day 20,200 (about 55 years). According to 
Munk, the streamlines should represent an exponentially 
decaying series of vorticies centered on the latitudinal 
axis. The vorticies are present across the width of the 
domain but do not appear to decay as expected. Truncation 
error is a probable cause for insufficient decay in this 
area. In the gyre north of 45N, the maximum value of the 
vertically averaged meridional current was on the order of 
10 cm sec In the southern gyre the maximum value of the 

same quantity was about 15 cm sec The magnitude of the 

total (vertical shear plus vertical mean parts) meridional 
currents in the northern and southern gyres were 10 and 35 
cm sec respectively. The higher latitudes in the model 
were more barotropic and the vertical mean current domina- 
ted, while the vertical shear current dominated in the 
highly baroclinic areas of the central portion. In the 
north the vertical shear current was about zero, while in 
the south the vertical shear current was about 20 cm sec 
directed northward , 

Figure 22 shows the sea surface temperature field 
at day 20,200. The mean meridional temperature gradient 
over the latitudinal domain was 21,4 degrees per 6000 km. 
The temperatures at both the north and south boundaries had 
decreased from their initial values. Upwelling along the 
central portion of the northern boundary was responsible 
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Figure 21. Ocean stream function obtained at 
from 400 km grid combined model. 



day 20,200 
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for this cooling at the surface. The small scale spatial 
fluctuation in the temperature field in low latitudes ap- 
pears to be a result of the linear viscous solution in 
which there is alternating northerly and southerly flow and 
therefore northerly and southerly advection of heat. 

Figure 23 shov/s the zonal winds at three levels in 
the atmosphere. The curve labeled (a) depicts the zonal 
mean surface wind. The easterly, westerly, easterly pattern 
from north to south is qualitatively like the patterns in 
the real atmosphere. However, the latitudes at which the 
direction of flow change occur are slightly north of the 
actual latitude of change. Curve (b) is the zonal wind at 
750mb and curve (c) is that quantity at the 250 mb level. 

The increase of westerlies with height in the central 
latitudes is in thermal wind balance with a temperature 
structure of warm air in the south and cold air in the 
north , 

While the structure of the atmosphere-ocean system 
that evolved from integrating the model equations over the 
400 km grid was not extremely accurate, it did suffice for 
initializing the integration of the fine grid model. 

2 . The 200 km Grid Combined Model Experiment 

The fine grid or 200 km grid part of phase I was 
initialized with data from day 20,200 of the 400 km grid 
model after spatial interpolation. The model equations 
were then integrated in time for 8700 days (approximately 
24 years). The adjustment to the fine grid was evident in 



108 



Laf if ude 



72 



18U— JL 

-40 



Figure 23. 
grid model 




Zonal atmospheric winds in m sec from 400 km 
at a) surface; b) 750 mb; c) 250 mb. 
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a slight initial decrease of the oceans volume averaged tem- 
perature. After a few hundred days, this temperature start- 
ed increasing and at day 21,000 the sea was warming at a 
rate of 0.8 degrees per century. This warming rate was main- 
tained for over 1000 days and then the rate decreased grad- 
ually to 0.5 degrees per century at day 28,900, the last day 
of the experiment . This warming rate was caused by a down- 
ward heat flux across the air-ocean interface of about 2.3 
-2 -1 

cal cm day . While the experiment could have been con- 
tinued until no heat flux was observed, it was believed that 
the large amount of computing time required to accomplish 
this would not have been worth the slight change in the 
final results. 

At the end of the 400 km grid model experiment, the 

total disturbance energy in the atmosphere was statistical- 

2 -2 

ly steady, oscillating about 50 m sec . During the first 

few hundred days of the fine grid integrations the average 

2 -2 

disturbance energy jumped to about 70 m sec , but then 
quickly returned to the statistically steady state which 
had previously been observed. The initial cooling of the 
ocean imparted an added amount of heat to the atmosphere 
in the southern half of the region generating more zonal 
potential energy which the baroclinic processes converted 
into disturbance energy. After the cooling ended, the 
warming of the ocean apparently occurred rather uniformly 
in space resulting in the generation of about the s arae 
amount of available potential energy as had been generated 
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prior to reducing the grid size. More discussion of the 
energy balance in the fine grid model is in a later section. 

The mean state of the atmosphere over the last 100 
days of this phase of the experiment is shown in Figures 24 
to 26. These maps along with Figures 28 to 30 that will be 
discussed shortly are deviation values of the height fields 
at various levels obtained from the geostrophic condition. 
The zero contour and plus and minus deviation contours are 
drawn. These maps are displayed in this way because the 
deviation of the stream function is predicted by the quasi- 
geostrophic atmospheric equations. The quantity appears 
as in all terras, and can be represented as V (ij; ’ (x ,y , t ) 

+ ^(t)), a deviation and mean part of the stream function. 
Therefore, any mean ip could be added to (or subtracted 
from) the deviation quantity. 

Figure 24 shows the extrapolated distribution of 
the 1000 mb height deviation analyzed at 12-meter intervals. 
There is a semi-permanent mid-latitude cyclone over the 
eastern ocean and a high at the same longitude south of 
the low. A high is present in the northernmost section of 
the region. This is worthy of note since there is no cold 
land mass modeled into the system to simulate a "Greenland” 
type effect. Figure 25 shows the distribution of the 250 
mb height deviation analyzed at 60-meter intervals. The 
flow is zonal throughout the central portion of the domain 
with a cyclonic polar vortex present in the north and a 
weak ridge present at the same longitude in the south. 
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Figure 25, Distribution of 100 
deviation in phase I. 
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Figure 26 shows the distribution of the 750 mb height de- 
viation analyzed at 30-meter intervals* The flov; is again 
zonal throughout the central portion of the region. How- 
ever, the lowest height is now in the northeastern corner 
of the area and the low latitude high is present in only 
the southernmost part of the domain. 

The vertical structure of the standing disturbance 
is such that the 1000 mb high in the polar regions lies 
belov7 the 250 mb low the same as in the real atmosphere. 

The warm low pressure area at the southern boundary at the 
1000 mb level becomes an area of high pressure at the 250 
mb level . 

The zonal wind averaged over the s ame 100 day per- 
iod is shown in Figure 27. The line labeled (c) shows a 
westerly flow at all latitudes at the 250 mb level with the 
jet located at 45N, coinciding with the strong height 
gradient in that region as seen in Figure 24. Curve (b) 
shows the weak polar easterlies at 750 mb level to the 
north of the low at 60N. The remainder of the flov7 is 
westerly with a weak jet at about 45N. Line (a) shows the 
zonal flow at the 1000 mb level. The polar easterlies 
extend from the northern boundary to a latitude just north 
of the cyclonic activity. Westerlies extend from that 
point to the latitude of the 1000 mb highs and easterlies 
exist to the southern boundary. A comparison of Figure 27 
with Figure 23 shows that the mean zonal winds did not 
change significantly after reducing the grid size. This 
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Figure 26 . Distribution of 100 day mean 
deviation in phase I. 
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indicates that for the atmospheric model, a 400 km grid 
would have sufficed with little loss of resolution. 

Figure 28 shows the distribution of the 1000 mb 
height contours at 25 meter intervals at day 28,900. A low 
dominates the western half of the region at this time and a 
high covers the eastern half. This figure was compared to 




Figure 27. Zonal atmospheric winds in m sec~^ from 200 km 
grid model at a) 1000 mb; b) 750 mb; c) 250 mb. 

the series of 1000 mb figures presented by Phillips (1956) 
and the gross features are almost identical to the early 
stages of his model atmospheric development. The main sur- 
face troughs and ridges are orientated so as to lean back 
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toward the northwest in the northern half and slightly to- 
ward the southwest in the southern half of the region. The 
convergence of momentum in the central portion increases 
the westerlies and accounts for the mid-latitude jet ob- 
served in this area (Figure 29). 

Figure 29 shows the distribution of the 250 mb 

height contours at 60 meter intervals for day 28,900 (79 

years). A wave of moderate amplitude is superimposed on 

the zonal flow with a trough present over the western half 

of the region and a ridge dominating the eastern half. The 

polar vortex observed in Figure 25 exists at this time 

also. The single wave structure of the atmospheric model 

forces a pattern similar to this to be found at any given 

instant. The disturbance moves and its amplitude varies 

with time from a zonal pattern to a highly amplified wave 

pattern. At this time the disturbance kinetic energy had a 
2 -2 

value of 54 m sec , This is quite large compared to the 
disturbance kinetic energy of the more zonal flow in Figure 
25. However, at other times in the experiment higher 
values of this quantity were observed. It is believed that 
the atmospheric state shown in Figure 28 is representative 
of a flow in which a wave disturbance of average amplitude 
is imbedded. At times the amplitude of the wave v/ill be 
greater than it is at day 28,900 and at other times less. 
Also, the wave will propagate across the domain and yield 
a time mean upper level flow as shown in Figure 25. 
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Figure 28 , Distribution of 1000 mb height deviation at day 
28,900 in phase I. 
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Figure 29. Distribution of 250 mb height deviation at day 
28,900 in phase I. 
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Figure 30, the distribution of the 750 mb height 
contour at 30-iueter intervals at day 28,900 , shows a closed 
cyclonic circulation directly beneath the trough at 250 mb. 
The amplitude of the wave motion is about the same as at 
250 mb and a closed high is present beneath the upper level 
ridge • 

The suggestion that the single wave in the x-direc~ 
tion would produce a series of cyclonic scale disturbances 
moving over the d oma in of the atmosphere-ocean model was 
born out by these figures plus examination of preceding 
maps. The wave moves eastward with a speed of about 600 km 
per day. This is considerably slower than the wave motion 
in Phillips’ experiments. 

The vertical heat flux was discussed earlier as be- 
ing proportional to the atmosphere-ocean temperature differ- 
ence. However, due to the method by which the sensible 
heat exchange is calculated, the vertical heat flux is 
actually proportional to the difference between the atmos- 
pheric surface temperature and the Fourier analyzed sea 
surface temperature. This temperature difference, averaged 
over the last 100 days of the experiment, is shown in 
Figure 31. In general, there is cooling of the atmosphere 
in the western sector and warming of the atmosphere in the 
east. While this heating may be representative of what 
occurs in the northernmost latitudes of the real atmosphere, 
it does not accurately portray what occurs throughout the 
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Figure 30. Distribution of 750 mb height deviation at day 
28,900 in phas e I . 
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Figure 31, 100 day average difference between the Fourier 

transformed sea surface temperature and the atmosphere sur- 
face temperature computed in phase I. Heating of the 
atmosphere occurs where values are positive. 
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majority of the oceanic regions where cold atmospheric air 
masses move from over a cold continent to a position over 
warmer water and are heated from below. The lack of colder 
continents in the model and the use of cyclic boundary con- 
ditions makes it impossible to simulate this phenomena in 
this simple model. The use of the Fourier analyzed sea 
surface temperature also contributes to the unrealistic 
vertical heat flux pattern of the model. The strong tempera- 
ture gradient observed in the real ocean along the east 
coast of continents and observed in the model along the 
western boundary of the ocean in Figure 34, had little 
effect on the heat flux pattern that evolved. The Fourier 
analysis smoothed this gradient and although the gradient 
did influence the new sea surface temperature somev/hat, it 
did not have the same influence observed in the real atmos- 
phere-ocean system. The Fourier analyzed sea surface tem- 
perature was used in the heat exchange to be consistent 
with the simple one-wave formulation of the atmosphere. 

The cyclic boundary condition forces the air to pass from 
a region of warm water on the eastern boundary immediately 
to the region of cooler water in the western area. The 
overall effect that the resulting heat flux pattern had on 
the results of this experiment are hard to determine. The 
strong heating observed in nature on the western side of the 
ocean in the mid-latitudes does not occur. The overall 
effect may be to smooth out the effect of the heating on 
the atmospheric temperature disturbance (only disturbance 



123 












\ M 











quantities are of interest since domain mean temperatures 
of the atmosphere do not change) and cause a more zonal 
pattern. Figure 32 shows the time mean atmospheric surface 
temperature averaged over the last 100 days of the combined 
model experiment. The temperature in this figure has a 
more zonal appearance than was expected after examining the 
mean 1000 mb height deviations (Figure 24). 

Figure 33 is the 100 day time mean of the stream 
function governing the vertical mean currents. The sub- 
arctic gyre, visible at about 55N in this figure, produces 
southward flow along the western boundary at that latitude 
and northward return flow 4-600 km east of the boundary. 

The subtropical circulation at about 35N has northward flow 
at the western boundary and southward flow 4-600 km to the 
east. The convergent flow in the boundary causes a sharp 
meridional temperature gradient at about 43N and contri- 
butes to the low temperature near 48N in the western bound- 
ary . 

Figures 34, 35, and 36 show the 100 day time means 

of the ocean temperature at the surface, 500 meters and 
1300 meters, respectively. These maps are vertically con- 
sistent and the low temperatures just described near 48N can 
be seen at all three levels. The thermal ridge observed 
in the northwestern sector of Figure 34, and also in Figure 
35, is caused by the nor thwes tvrar d flow on the east side 
of the subarctic gyre and somewhat by the upwelling along 
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Figure 35. Temperature pattern at a depth of 500 meters in 
ocean in phase I (C) , 
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Figure 36. Temperature pattern at depth of 1300 meters in 
ocean in phase 1(C) . 
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the western boundary seen in Figure 37, which contains the 
time averaged vertical motion Field at the depth of 300 
meters. In the fine grid, as v;ell as in the coarse grid 
discussed earlier, the flow in the western boundary regions 
in high latitudes is composed almost entirely of the verti- 
cal mean current. Here the flow is of this character and 
moves across the isotherms, advecting warmer water toward 
the northwest. The large area with horizontally isothermal 
water of about 5C located west of the ridging in the north- 
west quadrant of the ocean domain can be compared qualita- 
tively to a similar feature found east of Newfoundland in 
the North Atlantic, Sverdrup et al (1942) presented sea 
surface temperature maps that contained this feature in 
both winter and summer. 

In the northeastern quadrant the flow is weak and 
closely parallels the isotherms. Here, in the absence of 
appreciable vertical mean current, the vertical shear cur- 
rent is dominant and when the flow reaches the northeastern 
boundaries large sinking motions occur (see Figure 37) pro- 
ducing warm temperature at all depths there. Figures 34, 

35 and 36 all show warmer water in the northeast corner 
due to this phenomenon. The isotherm pattern produced by 
the model in this area can be likened to the real sea sur- 
face temperature pattern off the west coast of the United 
States. The area of warmest surface water in Figure 34 is 
observed in the southwestern corner of the domain. This 
area is analogous to the Sargasso Sea with its warm tem- 
peratures and weak currents. 
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The temperatures along the western boundary are fur- 
ther affected by the upwelling that takes place due to the 
thermal circulation. The central latitude is a region of 
weak eastward flow in the surface layer and westward flow 
at the lower layers . The moderate upwelling brings colder 
water to the surface, lowering the surface temperature 
along the boundary. An area of weak sinking motion in the 
south central portion of the domain in Figure 37 brings 
about the warm 8C water at the 500 meter level (Figure 35). 

From the stream function field in Figure 33, the 
vertical mean current can be computed. The subarctic gyre 
has a vertical mean current on the order of 20 cm sec 
The subtropical gyre has a vertical mean current on the 
order of 15 cm sec Table 3 shows the total meridional 

velocity, an estimate of the vertical mean meridional cur- 
rent as calculated from the stream function and the verti- 
cal shear meridional current calculated as a residual, in 
the western boundary current. The vertical mean current 
is almost entirely responsible for the meridional current 
in the boundary layer associated with the subarctic gyre, 
while the vertical shear current is dominant in the sub- 
tropical gyre. The maximum value of 51 cm sec ^ in the 
boundary current, representing a Gulf Stream type flow, is 
less than that found in the real world oceans where 100 
cm sec ^ have been observed for this type of a boundary 
current . 
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Latitude 


Total Merid. 
Current 


Vertical 

Mean 


Vertical 

Shear 


71.1 


2 . 6 


2 . 8 


-0.2 


67.5 


12.8 


5.6 


7.2 


63.9 


9.9 


10.0 


-0.1 


60.3 


-0.4 


-1 . 7 


1.3 


56 . 7 


-12.0 


-14 . 0 


2 . 0 


53.1 


-20 . 5 


-19.0 


-1.5 


49 . 5 


-18 . 7 


-14.0 


-4 . 7 


45.9 


-5 . 6 


-6.5 


0.9 


42.3 


11.1 


0.8 


10.3 


40.5 


17 . 7 


3.8 


13.9 


36.9 


37.5 


12 .5 


25 . 0 


33.3 


51 . 1 


15 . 5 


35 . 6 


29 . 7 


41.0 


8 . 0 


33.0 


26 .1 


20.5 


2.5 


18.0 


22.5 


6 . 7 


-5.1 


11 . 8 


18.9 


1.9 


-0.7 


2.6 



TABLE 3. Total meridional current (cm sec ) as produced 
by the combined model; vertical mean current as computed 
from the model stream function; vertical shear current 
calculated as a residual. 
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Figure 38. 100 day mean stress pattern produced by com- 

bined model . 

The stress acting as a driving force on the ocean 
was continuously varying as the tir.ie integrations progressed. 
A 100 day time mean zonal stress was calculated at a time 
six years prior to the end of the integrations and is shown 
in Figure 38. This stress was used as a constant force 

acting on the ocean during the ocean-only phase of the 

experiment. Now it can be interpreted as the mean stress 

responsible for the stream function field observed in 
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Figure 33. The curl of the stress has maxiina at about 3AN 

and SON which are the approximate positions of the centers 

of the subtropical and subarctic gy r e s , r e s pe c t ive ly . The 

positions of the maxiina are about 5 degrees further north 

than actually observed in nature. This is responsible for 

the northward shift of the gyres and their associated 

boundary currents. The stress varies from 1.20 to -1.15 
-2 

dynes cm . These are reasonable values lending justifica- 
tion to the stress formulation discussed earlier. 

The state to which the combined model progressed 
was acceptable in its resemblance to the real atmosphere- 
ocean system. While it was not the intention of this work 
to develop an operational model for forecasting either 
atmospheric or oceanic parameters, it is comforting to see 
a reasonable distribution of the various parameters. 

3 . The Ocean-Only Experiment (Phase II ) 

The zonal mean stress (Figure 38) that had been 
averaged over the 100 day period from day 26,650 to 26,750 
was used as a constant forcing function for the subsequent 
six year period over which the ocean-only model equations 
were integrated. The atmospheric surface temperature, T^, 
was also averaged over this same 100 day period. Figure 
39 shows a plot of the zonal mean temperature as a function 
of latitude. The temperature varied from 2.7 degrees 
centigrade in the north to 23.3 degrees in the soutli. This 
temperature was used as an equivalent atmospheric tempera- 
ture to which the sea surface temperature was compared to 
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determine the vertical heat flux. Except for the fixed 
atmospheric surface temperature which was independent of x, 
the same formulation was used for the heating as v/as used 
in the combined model. Now, however, one degree of freedom 
was removed and the ocean in a very short time became 
steady, except for a slow secular change. 




Figure 39. 100 day time mean atmospheric temperature (C) 

as a function of latitude. 



The instantaneous fields of the ocean dependent 
variables from day 26,700 were used as initial conditions 
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for the ocean--only phase of the experiment. During the six- 
year ocean-only experiment, data fields for all oceanic 
variables vere saved every 100 days and selected daily para- 
meters were saved at the same five points as previously 
discussed. 

The 2200 day (six year) integration of the ocean- 
only equations proceeded at a rapid pace. No problems were 
encountered with adjustment to the new constant forcing 
functions and a quasi-steady state was achieved. At the 
start of the 2200 day phase of the experiment the volume 
average temperature of the ocean was increasing at about 
0.5 degrees per century as previously discussed. Through- 
out the six year period the v;arming rate remained about the 
same, with a 0.5 degree per century v;arming still present 
at the end. This makes this author believe that perhaps 
this is as close to an equilibrium state as could be 
reached with the model as formulated, without the expendi- 
ture of an excessive amount of computing time. An examina- 
tion of the sea surface temperature. Figure 40, and the 
area average of the temperature at that level shows that 
the sea surface temperature is not changing at all. The 
rate of change of average temperature at the upper level 
is exactly zero. However, the average temperature at the 
lowest level is changing at a rate of 0.6 degrees per 
century, showing that the diffusive mechanism by which the 
deeper water must be changed is very slow. 
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Latitude 




Figure 40. Sea surface temperature from ocean-only model 
at day 28,900. 
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Figure 40 shows the sea surface isotherms at day 
28,900. An examination of this temperature field at day 
28,800 was made and the two maps were identical. Therefore, 
it was considered unnecessary to derive t ime mean map s of 
this quantity. The pattern is similar to the sea surface 
temperature field at the end of phase I (Figure 34) in that 
a strong temperature gradient is present in the western 
boundary layer at about 35N. The isotherms throughout the 
central portion of the domain are slightly more zonal than 
in phase I where there is an x dependence in the forcing, 
but the meridional gradient in this area remained about the 
same. The warm core tropical water (above 23C) does not 
extend into the boundary layer but is visible as a narrow 
band of warm water just north of the southern boundary. 

The major difference in the temperature fields be- 
tween the two phases is in the northern half of the domain. 
The thermal ridge present in the northwest quadrant in 
phase I is observed only as a small perturbation in the 6C 
isotherm in this phase. The thermal trough to the west of 
this feature is similar to that present in the combined 
model with a region of water with temperature between 4 and 
6 degrees extending over a thousand kilometers. The remain- 
der of the isotherms in the north have a more zonal appear- 
ance. There is a small ppcket of water below 2C in the 
northernmost region in phase II that is not present in 
phase I. The isotherms in the northeast sector vary signi- 
ficantly between the two phases. In pliase I the thermal 
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current flows into the northeast corner v/hen the strongest 
sinking motion occurs. . In phase II the vertical shear cur- 
rent flows into the eastern wall on a more eastv/ard heading. 
The maximum sinking motion is located about 800 to 1000 km 
south of its position in phase I. This can be verified by 
comparing Figure 41, the vertical velocity at a depth of 
300 meters in the ocean-only case, with Figure 37, the same 
quantity in the combined model. The northward curvature of 
the isotherms in the central portions of the eastern bound- 
ary is present in both phases, but is less pronounced in 
the ocean-only case. An examination of the meridional cur- 
rent in this area showed that there was weak northward 
motion in the northeastern boundary region in both phases 
which when coupled with the downwelling in that area, was 
responsible for the warm water there. The temperature 
patterns for the levels at 500 and 1300 meters resemble 
their counterparts from phase I (Figure 35 and 36) with the 
same basic variations as were seen in the upper level. The 
maximum in the temperature in the northeast corner has moved 
south about 800km as the region of sinking motion moved. 

The isotherms show no unexpected results and will not be 
depicted here . 

The vertical motion patterns play an important part 
in the adjustment of temperature patterns in the boundary 
regions. Strong vertical motions occur there due to the 
thermal circulation. A region of sinking water near AON 
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Figure 41. Vertical motion in cm day~^ at 300 meters at 
day 28,900 of phase II. Positive indicates upward motion. 
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in the western boundary region causes the weak gradient be- 
tween the 12C and 14C isotherms. The isotherms along the 
western boundary are shifted southward. That is, the 14C 
isotherm in the ocean-only model is at the latitude at 
which the intermediate 15C line (not drawn) was in the com- 
bined model . 

Figure 42 shows the stream function field for the 
ocean-only model at day 28,900. It was also found unneces- 
sary to compute a time mean of this quantity. The subarctic 
gyre and the subtropical gyre both appear to be slightly 
further north in the ocean-only model than they are in the 
combined model (Figure 33) • Those gyres have about the 
same intensity in the western boundary region as the corres- 
ponding gyres in phase I. This was expected since a mean 
stress from the combined model was their driving mechanism. 
The northwestward flow on the east side of the subarctic 
gyre is still present and causes some warm water to be 
advected, but at this time the large amplitude ridge that 
was observed in phase I was not present. The ocean-only 
sea surface temperature appears to be strongly tied to the 
imposed zonal atmospheric temperature, except in the west- 
ern boundary region and to a lesser degree in the eastern 
boundary region. 

The maximum value of the meridional current is re- 
duced somewhat in the ocean-only model. In the western 
boundary region, the maximum meridional velocity is 47 cm 
sec ^ in the ocean-only model compared to 51 cm sec ^ in 
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Figure 42. Stream function pattern at day 28,900 in 
phase II. 
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the combined model. Since the vortical mean current appears 
about the same in both models, the difference is in the 
vertical shear current. The difference in the predicted 
ocean variables between the ocean-only and combined model 
phases is small. Radical changes did not occur. The dif- 
ferences that are present are due to the removal of the 
time varying atmospheric forcing mechanism in both the ex- 
change of momentum and heat in the ocean-only model. 

4 . The A t mo sp he re -Only Model (Phase III ) 

. ■ y 

The zonal part, T , of the combined model sea- 

sea 

surface temperature was averaged over the same 100 day 
period as were the zonal stress and atmospheric surface 
temperature. This quantity was used as the underlying tem- 
perature with which the atmospheric temperature was com- 
pared for the calculation of the vertical flux of sensible 

heat. Figure 43 shows the profile of the time and longi- 

1 

tudinal mean of sea surface temperature, T , which 

sea 

closely resembles the mean atmospheric temperature shown 
in Figure 39 . 

The initial conditions for phase III were the 
Fourier coefficients, A through F, from day 26,700 of 
phase I. There was essentially no adjustment required and 
the 2200 day integration was conducted, with data retained 
as before, in a minimum of computing time. At the end of 
this phase the mean state of the atmosphere over the last 
100 days was examined along with the state of the atmosphere 
from day 28,900. These maps were compared with the atmos- 
phere at the end of phase I. 




Figure 43. 100 day mean zonal component of sea surface 

temperature (C) . 

Figure 44 shows the 100 day mean distribution of the 
1000 mb height deviation as a function of latitude. Since 
the thermal forcing is now a function of y only, the mean 
fields are purely zonal (with the exception of small errors 
due to averaging and truncation). The deviation in height 
extends from 18 meters at 35N to -28 meters at 54N and the 
deviation pattern is such that a geostrophic wind regime 
of easterlies, westerlies, easterlies would exist (as it 
does). Figure 45 shows the 100 day mean distribution of 
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Figure 44. Distribution of 100 day mean 1000 rab height 
deviation as a function of latitude. 

the 750 rab height deviation as a function of latitude. The 
deviation goes from 118 meters at the southern boundary to 
-120 meters at the northern boundary. The tight negative 
height gradient through the central latitudes agrees with 
the existence of a regime of westerlies in that area. 
Figure 46 depicts the distribution of the 100 day mean 
distribution of the 250 mb height deviation as a function 
of latitude. The gradient that exists from 410 meters in 
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the south to-407 meters in the north also agrees with the 
westerly regime at that level. 




Figure 45. Distribution of 100 day mean 750 mb height 
deviation as a function of latitude. 

Instantaneous maps at day 28,900 of the distribution 
of the height deviation show a highly disturbed atmosphere. 
Figure 47, the 1000 mb height deviation field at day 28,900 
analyzed at 25-meter intervals, shows a pattern that has 
a disturbance amplitude which is about the same as its 
counterpart in the combined model (Figure 28). Figure 48 
shov7s the distribution of the 250 mb height deviation 
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analyzed at 60 rneter intervals. This flow is quite similar 
to that at the same level at the same time from the combined 
model (Figure 29). Figure 49 depicts the distribution of 
the 750 mb height deviation analyzed at 30 meter intervals. 
The wave-like structure here has a greater amplitude and 
the deviation is larger than the same fields in the combined 
model (Figure 30). The occurrence of the difference in the 
amplitude of the patterns is not surprising since the model 
atmosphere is continuously changing. It v/ould be a coin- 
cidence if the structure of both the atmosphere-only and 
combined model atmosphere were identical at any given time. 




Figure 46. Distribution of 100 day mean 250 mb height 
deviation as a function of latitude. 
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Figure 47. Distribution of 1000 mb height deviation at 
day 28,900 in phase III. 
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Figure 48. Distribution of 250 mb height deviation at 
day 28,900 in phase III. 
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Figure 49. Distribution of 750 mb height deviation at 
day 28,900 in phase III. 
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The 100 day average zonal wind profile in the atraos- 
phere-'only model is very similar to that in the combined 
model shown in Figure 27, A very slight reduction in the 

jet maximum at 250 mb is the only difference so the atraos- 

phere-only wind profile will not be presented graphically. 

The time averaged results of the atmosphere-only 
phase of the experiment were very similar to the time 
averaged atmospheric state arrived at by the combined model. 
Since the sea surface temperature is a slowly varying para- 
meter compared to the atmospheric time scale, holding it 

constant approximates long term actual conditions fairly 
well. It would be surprising if a radical change occurred 
in the time averaged fields in the atmosphere-only formu- 
lation. 

5 . Energy balance 

An examination of the energy in the atmosphere was 
made at the completion of phase I and phase 111. Rather 
than study instantaneous values of the energies and energy 
transformation as was done in the earlier discussion of the 
energy balance, values averaged over the last 100 days of 
the experiment were examined. The energy flow convention 
as established in Figure 4 was used. 

The mean energy exchanges for the combined model 
and the at mosphere-'only model are shown in Figures 50 and 
51, respectively. It is interesting to note the similarity 
in the flow diagrams in the two phases. The diabatic heat- 
ing generates about the same amount of zonal potential 
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Figure 50. Flow diagram of mean atmospheric energy from 
the combined model. Energy and energy transformation 
units are in watt m"^ ( = 10^ erg cm“ sec”^). 

-2 

energy in both cases, 1.16 watt m in phase I and 1.19 

watt m in phase III. The conversion of zonal potential 

to eddy potential and then into eddy kinetic energy is 

-2 

identical in both phases at a value of 1.26 watt m . The 
conversion of eddy kinetic to zonal kinetic energy is also 
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Figure 51. Flow diagram of mean atmospheric energy from 
the atmospheric-only model. Same units as in Figure 50, 

very similar in the two phases v;ith values of 0.37 and 0,33 
watt m , The dissipation of the kinetic energy (both 
zonal and eddy) by both internal and surface friction are 
also very similar. The only conversion that differs by a 
significant percentage is the conversion of zonal kinetic 
to zonal potential energy. However, this has been pointed 
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out earlier as being a very small quantity and no comment 
on the difference between the two phases is offered here. 

It is not surprising that these flow diagrams re- 
semble each other. The only difference in the two models 
lies in the method by wh ich heat is exchanged be tv/een the 
atmosphere and ocean or underlying wet surface as the case 
may be. In the calculation of the sensible heat flux 
necessary for the generation of zonal potential energy in 
the combined model only the zonal part of the sea surface 
temperature was used. Therefore, when the disturbance 
part of the sea surface temperature was neglected in the 
atmosphere-only model, it did not affect the generation of 
energy very much. In the other energy exchanges the differ- 
ences in the disturbance quantities were not sufficient to 
make much difference. 

These diagrams compare favorably with Oort’s dia- 
gram (Figure 6) as did the earlier energy diagram (Figure 
5). The dissipation of eddy kinetic energy is still too low 
and no new explanation is offered for this. It is somewhat 
surprising that the instantaneous values of the energy ex- 
change shown in Figure 5 compare so favorably with Figures 
50 and 51 in so many ways. The fact that the time rate of 
change of kinetic and potential energy (both zonal and 
eddy) were other than zero in that case shows that either 
there were small errors in calculations, or an exact steady 
state in the energies did not exist, as was stated in the 
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earlier discussion. This is in accord with present know- 
ledge of the real atmosphere which is never in exact steady 
state . 

B. SPECTRUM ANALYSIS 

This section will briefly discuss the spectrum analysis 
techniques used in this study. A computer program in the 
BIOMEDICAL series (BIMED 02T) developed at the UCLA Medical 
Center (Dixon 1970) was modified and utilized for all 
spectral analyses conducted on the data generated in this 
study . 

1 . Basic Definitions 

Prior to discussing the procedures used or the 
results obtained, a brief review of the definition of some 
of the basic statistical parameters is given. 

a. Variance - the mean square deviation from the 
mean of a variate. 

b* Standard deviation - the positive square root 
of the variance, used to measure the degree of variability 
of a variate. 

c. Auto-correlation coefficients - linear correla- 
tion coefficients between a time series and the same time 
series an interval of time later. 

d. Auto-covariance ^ an auto-correlation function 
that has not been normalized by the variance. 

e. Cross-correlation - correlation between two 
time series, one of which is at an interval of time later. 
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f. Cross-covariance - a cross-correlation function 
that has not been normalized by the variance. 

g. Power spectrum of time series - a function of 
frequency which gives the contribution of oscillations at 
each frequency to the variance of a time series. 

h. Cospectrum - a function of frequency which 
gives the contribution of oscillations at that frequency 
to the total cross-covariance of two time series. 

i. Quadrature spectrum - a function of frequency 
which gives the contribution of oscillations at that fre- 
quency to the total covariance of two time series X(t) and 
Y(t) obtained when all the harmonics of time series X were 
delayed by a quarter period but the Y series remained un- 
changed . 

j. Phase - the relative phase of the harmonics 
X(t) and Y(t) is measured by the ratio of the quadrature 
spectrum to the cospectrum. 

2 . Spectral Techniques 

Many investigators (for example, Yanai et al, 1968; 
and Wallace 1971) have utilized spectral techniques to 
study tropospheric phenomena primarily in areas of sparse 
data where a single station’s time varying record is 
extremely important. The same type of spectral techniques 
can be used to study the atmospheric and oceanic data ob- 
tained from the numerical integrations conducted in this 
study. At five selected points (Figure 20) six years of 
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data were recorded at 1.125 day intervals. At each point 
1960 observations were obtained. These data were analyzed 
and a power spectrum was obtained for most variables at all 
five points. Cross-spectrum analyses were conducted on sev- 
eral variables yielding cospectra, quadrature spectra and 
phase relationships. 

More attention was focused on point C because of 
its position in the western boundary region. It was located 
in an area where meridional velocities were strong, tempera- 
ture gradients were tight and upwelling was important. 
However, the other points in Figure 20 were not neglected 
and comparisons of the resulting spectral statistics between 
points were made. 

The data were recorded at an interval, (A of 

1.125 days (three oceanic time steps) during the three six- 
year experiments. That is, the data were available for 
analysis from the combined, ocean-only and atmosphere-only 
phases of the experiment. The decision to record data at 
1.125 days was a critical one. It would require an exces- 
sive amount of storage space to keep a longer record of 
data at shorter recording intervals and problems could be 
encountered if the interval was too long. It is important 
to ensure that the series yhich is being examined does not 
contain substantial amounts of variance from frequencies 
outside the range that cai'i be resolved with the spectra. 
Spectral information for frequencies less than one cycle/ 

2(A t)^, the Nyquist folding frequency, can be obtained from 
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the data. Frequencies higher than this are "folded” or 
"aliased" into the lower frequency region. Wallace (1971) 
points out that the high frequency limit does not appear to 
pose any serious problem from aliasing because there is 
little variance associated with frequencies higher than 0.5 
cycles per day in the intermediate to large scale disturb-- 
ances . From experimental work Wallace concluded that there 
were only minor differences in spectra computed v/ith daily 
or twice-daily data. Since the data examined in this study 
is of the daily variety, his conclusions, along with the 
fact that higher frequencies due to advection of short 
waves cannot exist in the model atmosphere leads the author 
to believe that the time interval chosen for this data re- 
cord was a sound one. 

In the BIMED 02T program, the power spectrum is ob- 
tained by performing a harmonic analysis on the auto-co- 
variance function of the quantity to be studied. After 
each spectral estimate was multiplied by the interval of 
frequency ( A f) of the band width, the spectral estimates 
for that frequency band can be summed to obtain the contri- 
bution of the fluctuations within that specific frequency 
band width to the variance of the time series which is 
under examination. 

3 . Results of Spectral Investigation 

The primary point of interest in the spectral in- 
vestigation was point C (Figure 20) in the western bound- 
ary region. A thorough investigation of atmospheric and 
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oceanic variables was conducted at this point with the fol- 
lowing results . 

a. Study of data at point C. 

The data obtained from phase I (combined model) 
of the experiment were examined first. The Fourier coeffi- 
cients from the atmospheric part of the combined model were 
analyzed and Figure 52 shows the power spectrum of A and B 
(solid line) and C and D (dashed line), obtained from Equa- 
tions (2.22) to (2.25), plotted on a linear scale against 
frequency (cycles per day) . The power spectra of A and B 
are identical as are the spectra of C and D. The amplitude 
of the power spectra of C and D is slightly less than 10% 
of that of A and B. There is a pronounced peak at about 
7.5 days in the power spectra for these quantities. This 
was expected since the phase speed of the imposed baroclinic 
disturbance in the combined model was estimated at 600 km 
per day, a phase speed that would yield a period of 8 days. 
In all likelihood, the estimate of the phase speed of the 
disturbance was slightly low. The 7.5 day period indicates 
that 640 km/ day was the speed with which the disturbance 
moved across the 4800 km wide domain. 

Figure 53 shows the power spectrum of atmos-' 
pheric disturbance kinetic energy. There are some small 
peaks in the spectra in the frequency band from .08 to .24 
cycles per day (12.5 to 4.2 days) but the major peaks are 
clearly at the lower frequencies. There are peaks at 
periods of about 32 and 62 days indicating cycles in the 
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Figure 52. Power spectra of A,B,C and D at point C in 
phase I. Period of spectral peak indicated in days. 
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Figure 53. Power spectrum of atmospheric disturbance 
kinetic energy in phase I. Period of spectral peaks 
indicated in days. 
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atmospheric energy on the order of one and two months. If 
seasonal variations had been included in the specification 
of the solar radiation and latent heat release, different 
frequencies of oscillation would probably be observed. The 
annual fluctuation in the above mentioned parameters would 
probably be reflected in the energy spectrum. 

Figure 54 shows the power pectrum of E, the zo- 
nal mean part of that is proportional to the mean height 

and related to the mean geostrophic wind in the atmosphere. 
The quantity has a large spectral peak at a period of 32 
days and smaller peaks at 25 and 62 days. An examination 
of the cospectrum and phase between E and the disturbance 
energy shows that the disturbance energy rise at the 32 
day period lags the increase in the spectrum of E at that 
period. This suggests that the mean flov; increased first 
and because of barotropic amplification transferred energy 
to the disturbance. 

Figure 55 depicts the power spectrum of the 
oceanic stream function. Peaks in the spectra are seen at 
periods of 7,5 days (small), and 32 and 62 days (large). 

This is the first quantitative indication of an air-sea 
interaction phenomenon in the model. The 7.5 day peak in 
the power spectrum of the stream function is caused by the 
contribution of the A and B coefficients to the disturbance 
component of the atmospheric wind stress (Figure 52). 

Since the disturbance components of the stress are small the 
fluctuations they cause are small and therefore the peak in 
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Figure 54. Power spectrum of E at point C in phase 
Period of spectral peaks indicated in days. 



J- 



0.40 



I. 



164 



42 




Frequency [cycles/day] 



Figure 55. Pov/er spectrum of ocean stream function at 
point C in phase 1. Period of spectral peaks indicated in 
days . 
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the spectruni is smaller at this frequency. The larger peaks 
at about 32 and 62 days are caused by the oscillations in 
the zonal component of the wind stress which is proportional 
to E. An examination of the cospectrum and phase between E 
and the ocean stream function shows that the increase in the 
mean atmospheric flow occurs 3.5 days prior to the fluctua- 
tion in the stream function at the 32 day period. Since E 
had spectral peaks at these frequencies it is not surprising 
that peaks in the same region of the spectrum appear in the 
ocean stream function. Figure 56 shows a 100 day time 
series of E and the oceanic stream function obtained by 
averaging the recorded data over seven day intervals. The 
averaging was done to remove the 7,5 day oscillation from 
the data to examine the longer period fluctuations in more 
detail. It is easy to see that the long period oscillations 
are present in both E and the stream function and that the 
maxima in E leads the maxima in the stream function by 
several days . 

Figure 57 shows the power spectrum of the zonal 

mean of the atmospheric surface temperature, T^ . Since T^ 

is proportional to F, the peaks in the spectrum are also 

representative of peaks in the spectrum of F. This was 

verified by spectrally analyzing F separately, the plot of 

which is not presented. There are peaks in T^ of almost 

equal amplitude at 7.5 and 32 days, and smaller peaks at 

periods of 25, 62 and 111 days. Figure 58 shows the power 

spectrum of T , the sea surface temperature. Small peaks 

sea 
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Figure 57. Power spectrum of atmospheric surface tempera- 
ture at point C in phase I. Period of spectral peaks 
indicated in days. 
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Figure 58, Pov;er spectrum of sea surface temperature at 
point C in phase I. Period of spectral peaks indicated 
in day s . 
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can be seen at periods of about 25 and 32 days and a larger 
peak at a period of about 62 days. These peaks coincide 
v/ith some of the peaks in the spectrum of , Figure 59 
shows the power spectrum of the vertical heat flux betv;een 
the atmosphere and the ocean. Peaks in the spectrum occur 
at 7,5, 32 and 62 days. The longer period peaks are at the 
same period as the peaks in the atmosphere and ocean surface 
temperatures, indicating an exchange of heat due to a tem- 
perature difference between the two media, A comparison of 
the time series of the atmospheric surface temperature, the 
vertical heat flux and the sea surface temperature obtained 
by averaging over seven day periods is shown in Figure 60, 

It is readily observed that oscillations on the order of 
30 days are present in all three series and not so obvious 
fluctuations on the order of 60 days are present in the sea 
surface temperature and the vertical heat flux. It can be 
seen that peaks in the sea surface temperature led peaks 
in the vertical heat flux which in turn led peaks in the 
atmospheric surface temperature. 

An examination of the cospectrum and phase be- 
tween E and T shows that the inci*ease in E at about 32 
sea 

days precedes the fluctuation in sea surface temperature by 

11 days. This V7as perhaps a result of an increase in the 

mean conal wind stress causing an increase in the wind 

driven part of the ocean circulation and an increase in the 

advection of heat from the south at point C. An examination 

of the cospectrura of T and the ocean stream function 

sea 
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Figure 59. Power spectrum of vertical heat flux at point C 
in phase I. Period of spectral peak indicated in days. 
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shows that the fluctuations in the stream function at 32 
days occurs 7.6 days prior to the oscillation in T 

sea 

These phase relationships coupled with the previously stated 
facts leads to the follov/ing conclusions. The oscillations 
which were reflected as peaks in the power spectra of the 
variables were caused by an air-sea interaction phenomenon 
and were not merely spurious reflections of energy due to 
sampling or calculating problems. An explanation of the 
probable series of events causing these interactions is as 
follows. The increase in the mean atmospheric flow, re- 
flected in E, caused an increase in the wind driven part 
of the ocean circulation 3.5 days later. This increased 
circulation brought about an increase in the advection of 
heat from the south at point C and an oscillation in the 
sea surface temperature, at a period of 32 days, occurred 
7.6 days later. The total time from the increase in E to 

the increase in T was about 11 days. The increase in 

sea 

T caused an increase in the vertical heat flux and an 

sea 

increase in the atmospheric surface temperature. An exam- 
ination of the cospectrum and phase of and T^ shows 

that the oscillation in T at 32 days leads the rise in 
the spectrum of T^ by 6.2 days. An exchange of energy be- 
tween T^(F) and E is then required to maintain the cycle. 

An examination of the cospectrum and phase of E and F shows 
that the peak in F occurs prior to the peak in E at 32 days. 
When these interactions are taking place xn the combined 
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model and oscillations at various frequencies are present, 
almost every dependent variable in the model will be affec- 
ted in some way. It has been shown that the long period 
oscillations which are present in the coupled model include 
the ocean as an important part in the oscillation cycle. 

It should not be surprising to find that these long period 
oscillations do not exist when either the atmosphere or 
ocean is held constant in time. 

The data from phase II of the experiment were 
examined next. With the time varying forcing functions 
removed, the ocean became steady in a short time with no 
significant peaks in any of the spectra. Figure 61 shov/s 
the power spectrum for This is typical of all the 

ocean variables in this phase in that the power in the 
spectrum is in the lowest frequencies. A comparison of 
Figure 61 and 58 provides a convincing argument that the 
time varying atmospheric circulation is responsible for the 
transient responses of the ocean surface temperature at 
periods of about one and two months. 

An interesting, and initially perplexing, peak 
appeared in the ocean stream function at 7.5 days in the 
ocean-only data (not shown). The ocean-only model required 
a period of time to adjust to the new forcing and during 
that time the stream function oscillated and then became 
constant. Although the energy in the spectrum was diluted 
by the long period of constancy, the initial oscillation made 
itself evident at the frequency at which it fluctuated. 
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Figure 61. Power spectruin of sea surface temperature at 
point C in phase II. 
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The data from phase III (atmospherC'-only ) were 
analyzed and Figure 62 shows the power spectra of A and B 
(solid line) and C and D (dashed line). The four quantities 
all have a peak at a period of 8.4 days. This is different 
than theperiod at which there is a peak in the spectra of 
the same quantities in the combined model. For the 4800 km 
wavelength used in this study, a period of 8.4 days yields 
a phase velocity of 570 km day This is about 12% slower 

than the phase speed for the atmospheric wave in the com- 
bined model. By eliminating the zonal variation in the sea 
surface temperature, the sensible heat exchange was changed. 
A comparison of the combined model and a tmos pher e -only model 
shows that the thermal wind was reduced in the latter. 

This reduction in the thermal wind, shown in Table 4, was 
influential in retarding the phase speed of the disturbance. 

The power spectrum of the disturbance kinetic 
energy is shown in Figure 63. A large peak at the 6.7 per- 
iod is observed along with a smaller, wider peak in the 45 
to 50 day period range. The peak at the 6.7 day period, 
which is near the 8.4 day period in the disturbance quanti- 
ties, is perplexing in light of the fact that there is no 
sharp peak around the 7.5 day peak in the disturbance kine- 
tic energy (Figure 53) in the combined model. A possible 
explanation is that the 6.7 day peak is very sharp and 
while it has a considerable amplitude, the total contribu- 
tion to the variance of the disturbance is not much greater 



176 



2800- 



2000 - 



“D 

a 

u 




u 

o 



(N 

E 

u 

1200 - 




400 



C D 




0.08 0.16 0J24 



032 



0.40 



Frequency [cycles/day] 

Figure 62. Pov/er spectrum of A,B,C and D at point C in 
phase III. Period of spectral peaks indicated in days. 
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LATITUDE 


THERMAL WIND 

(I) 


THERMAL WIND 
(III) 


72.0 


0.0 


0 . 0 


68.4 


1,2 


3.0 


64 . 8 


4.2 


3 . 7 


61 . 2 


3.9 


3.5 


57 . 6 


4.6 


4.3 


54.0 


5 . 7 


5.6 


50.4 


6.6 


7.5 


46 . 8 


8.5 


9.4 


43 . 2 


10.2 


9.2 


39.2 


9.1 


7.5 


36.0 


6.8 


5.6 


32.4 


4.4 


3.9 


28.8 


3.0 


3.1 


25 . 2 


2.3 


2 . 8 


21.6 


1.3 


2.5 


18.0 


0.0 


0.0 



TABLE 4. Comparison values of mean atmospheric thermal 
wind from combined model (I) and atmosphere-only model 
(III), (in m sec"l). 
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Figure 63. Power spectrum of atmospheric disturbance 
energy in phase HI. 



179 



than the small diffuse indication of a peak around the 7.5 
day period in the combined case. In the spectra for the 
quantities E and F, peaks at these same periods are ob- 
served, showing the exchange occurring between the disturb- 
ance and mean flow. 

An examination of data in phases I, II and III 
at point C indicates that the peaks in the spectra noted at 
about one and two months in the combined model depend on 
atmosphere-ocean interactions for their existence and are 
not simply spurious or caused by aliasing. The fact that 
these interactions are on the order of one and two months 
is not the most significant point, because .the period 
would undoubtedly change with some other period of atmos- 
pheric forcing. What is significant is that these longer 
period fluctuations can be attributed to air-sea interac- 
tions. The absence of such peaks in the spectra in either 
the atmosphere-only or ocean-only phases shows that the 
fluctuations causing the peaks are phenomena occurring only 
when the freedom found in the coupled model is allowed. 
Although there is a peak in both atmospheric and oceanic 
variables in the coupled model at 7.5 days, this is not 
considered as significant as the peaks at the longer per- 
iods because this frequency of oscillation was forced into 
the model with the selection of the 4800 km wavelength for 
the single wave in x. Changing this wavelength \ 7 ould cause 
the peak in the spectrum at 7.5 days to shift and undoubted- 
ly would have many other effects on the model. It would 
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be interesting to see if the peaks in the spectra at the 
lower frequencies would be shifted. 

Spectral analyses of the kinetic energy in the 
ocean of the combined model at the three ocean levels were 
made at point C. An approximation to the total kinetic 
energy at point C was used for this study. The total kine- 
tic energy should include the vertical mean kinetic energy 
and the vertical shear kinetic energy. A mathematical for- 
mulation for this is 

K-E = y[(u* - 1^)^ + (v’ + 1^)^] (6.2) 

ocean 2 9y 3x 

In recording the data at point C no provisions were made 

for the calculation of and so the zonal component of 

the vertical mean current at this point was neglected. The 

Zonal component of the vertical mean current was small and 

therefore this was not a bad assumption. The kinetic 

energy computed at each oceanic level was 

'^•’^ocean ■ 

Figure 64 shows the power spectrum of kinetic 
energy at 100 meters in the ocean of the combined model 
versus frequency in cycles per day, both plotted on log 
scales. Figure 65 is the same plot for 500 meters and 
Figure 66 shows the 1300 meter depth analysis. Peaks in 
the spectra at periods of 62, 25-32, and 7.5 days are ob- 
served at all three levels. This is consistent with earl-- 
ier findings. The amplitude of the 25-32 day period 
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Frequency Icpci) 

Figure 64. Power spectrum of oceanic kinetic energy at 100 
meters at point C in phase I. Periods of spectral peaks 
indicated in days. 
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Figure 65, Power spectriiin of ocean kinetic energy at 500 
meters at point C in phase I, Periods of spectral peaks 
indicated in days. 
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Frequency (cpcl) 

Figure 66. Power spectrum of ocean kinetic energy at 1300 
meters at point C in phase I. Periods of spectral peaks 
indicated in days. 
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is the largest at the upper level v/hile the amplitude of 
the 62 day period and 25-^32 day period are about the same 
at the two lower levels. A comparison of Figures 65 and 66 
shows that the spectra for the two lowest levels in the 
model nearly coincide. The maximum amplitude in the power 
spectra decreased by an order of magnitude between the 100 
meter level and the 500 meter level. This lessening in the 
variance of the kinetic energy between these levels is not 
surprising since the lower levels are less active. The 
amplitude remained about the same between the 500 and 1300 
meter levels . 

Thompson (1971) analyzed a long series of data 
recorded by the Woods Hole Oceanographic Institution at a 
site located at 39 20N, 70W in the North Atlantic. In the 

model this would be analogous to a position approximately 
200 km northeast of point C, The data were recorded at 
900 second intervals. The folding frequency associated 
with this small recording interval allowed Thompson to in- 
vestigate the shorter period oscillations. Some gaps in 
the data forced the estimation of some of the spectra in 
the low frequency region. He analyzed the northward and 
eastward components of velocity and through a summing and 
averaging process calculated the *^kinetic energy density^' 
at several depths in the ocean (not shown). None of these 
spectra showed any peaks at periods greater than one day, 
although os ciallat ions with periods of about ten hours were 
observed. The spectra for the three lowest depths were 
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nearly coincident, and there was a large decrease in the 
amplitude of the spectra between 100 meters and 500 meters. 
In this feature the study of Thompson and the study con- 
ducted in this v7ork were similar, 

b. Study of data at point A. 

At point A (northeast area) the power spectrum 

of the coefficients A,B,C and D were greatly reduced [a 

2 -12 

maximum value of 16 (cm sec ) cycles per day for A and B 
as compared to 675.0 at point C], although the peak in the 
spectra was at the same period, 7,5 days. An examination 
of the coefficient, E, shows the spectral peak at 32 and 62 
days and a large peak at 120 days. Looking back at Figure 
53, a smaller peak in the region of the 120 day period can 
be seen also. The oceanic stream function at point A shows 
oscillations that have smaller amplitudes in the power 
spectrum than do those of the same quantity at point C. 

The biggest difference between the power spectra of vari- 
ables at point A and those at point C occurred with T 

sea 

There ar^ no peaks at the 32 or 62 day periods and all the 
power is in the lowest frequency, much the same as seen in 
Figure 58, the same quantity in the ocean-only case. 

At point A in the ocean-only phase, the power 
spectra resembled those from the same phase at point C, 

The maximum of the power spectrum that existed in the low- 
est frequency was increased but will not be discussed fur- 
ther since it adds no new information. An increase in the 
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power in the lowest frequency nierely indicates that the long 

terra change v/as larger. In the a traospher C'-only phase, the 

power spectra for the Fourier coefficients resemble those 

at point C, except that the raaxiraura of the spectra at the 

8,4 day period is reduced to approximately 10% of the value 

it had at point C. At point A, the meridional temperature 

gradient is rather weak and the amplitude of the baroclinic 

wave is much less than in the lower latitudes of point C. 

c. Study of data at Point B, D and E 

At point B (central area) the power spectra of 

the coefficients A, B, C, and D were slightly larger than 

those at point C. This location at 45N is in an area of 

large thermal wind, so the amplitude should be greater. 

Except for the increase in amplitudes the spectra are very 

similar to those of the variables at point C (Figures 52 

through 56). The power spectrum of T , though, is more 

sea 

like that at point A, with all the power at the lowest fre- 
quency. This indicates extremely long period fluctuations, 
probably too long for this data record to portray. 

The ocean-only variables at point B resemble 
those at point A and will not be discussed. The Fourier 
coefficients, A through B , have a peak at the same 8.4 day 
period and the maximum in the power spectrum is increased. 

No new information was obtained from point B. 

A separate discussion of the data at point D 
V7ill not be made because of similarity of results to point 
B. Strange results were obtained when the data at point E 
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were spectrally analyzed. The Fourier coefficients, A 
through D, had multiple peaks in the spectra that V7ere not 
representative. It v/as decided that this point was selec- 
ted too close to the southern boundary and the artificial- 
ity of the wall caused unusual oscillations in the 
atmospheric data. The results from point E were discarded 
with no further comment. 

The results of the spectral study were that the 
ocean circulation caused transients to appear in the atmos- 
phereic circulation in the one and two month periods and 
possibly with a period of 120 days and that similar tran- 
sients appeared in the ocean circulation due to the pre- 
sence of the fluctuating atmosphere. A great deal of 
information has been wrested from the spectral statistics 
discussed here, but a. great deal more could be obtained if 
time permitted . 
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VII . SUMMARY 



A numerical model of a combined atmospheric-ocean system 
was developed for use in an investigation of the long-term, 
large-scale air-sea interactions and the effect of these in- 
teractions on transients that appear in the atmosphere and 
oceanic circulations. The problems associated with observa- 
tional studies of the air-sea system were eliminated by 
using a numerical model that could produce a long-term simu- 
lation of an interacting system in a reasonable amount of 
computing time. 

The atmospheric model used in this study was fashioned 
after the two level, quasi-geos tropliic model developed by 
Phillips (1956). New dependent variables were introduced 
by the use of a one-dimens ional (in x) Fourier series, in 
which a single mode was retained, to represent the origi- 
nal variable. This method of approach allowed the rapid 
integration of the atmospheric equations since the number 
of degrees of freedom were greatly reduced. Without this 
modification the combined model could not have been inte- 
grated for a sufficiently long time period to obtain mean- 
ingful results. 

The ocean model utilized in this work was a modified 
version of Haney^s (1971a) primitive equation model. The 
modifications included a reduction in the number of levels, 
a reduction in horizontal region size (and the number of 
grid points), a reduction in the depth of the basin and the 
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conversion to Q-plane (cartesian) geometry. The primary 
reason for these modifications was to increase the ratio of 
simulation time to computing time and permit a long term 
simulation which was accomplished. 

The coupling of the model was accomplished through the 
exchange of momentum through the wind stress and the ex- 
change of heat through the turbulent exchange process which 
was proportional to the air-sea temperature difference. 
Assumptions were made to simplify both of these exchanges. 
The atmospheric model, however, had to be restricted to 
simplified formulations if it was to be integrated success- 
fully for the necessary length of time. The absence of 
moisture in the atmospheric model required further simpli- 
fication in the formulation of the latent heat flux. These 
simplifications were indeed necessary and an examination of 
the results showed that the model produced atmospheric and 
oceanic circulations that were adequate for this experiment 
in their portrayal of the real atmosphere-ocean system. 

The experiment was conducted in three phases. Phase I 
was a long term integration of the combined model. The 
computing time was reduced by the use of a larger grid size 
for the first portion of this phase. The 400 km grid dis- 
tance was that used to accelerate the integration of the 
combined model equations in the first part of Phase I was 
just adequate to resolve the western boundary current. 

The coefficient of lateral diffusion of heat had been re- 
duced to prevent large horizontal transports of heat by 
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diffusion and keep the meridional temperature gradient at a 
realisticly high level. The lateral diffusion of momentum 
was reduced, but could not be made small enough so that a 
truly nonlinear solution could be obtained without having 
the width of the boundary current become smaller than the 
horizontal grid distance. The unrealistic spatial oscilla- 
tions seen in the results obtained with the 400 km grid were 
a direct result of the large grid size. 

The spatial oscillations observed in the 400 km grid 
model results quickly disappeared in the 200 km grid version 
of the combined model. A state of quasi-equilibrium was 
reached in the long terra integrations of the combined model. 
The results shown throughout this thesis were considered to 
be satisfactory for the purpose of investigating the effect 
of the oceanic (atmospheric) circulation on the time-depend- 
ent behavior of the atmosphere (ocean) . Some features in 
the results considered worthy of note are summarized below. 

A western boundary current developed with very strong 
northward flow, analogous to the Gulf Stream, on the west 
side of the sub-tropical gyre and strong southward flow, 
analogous to the Labrador Current, on the west side of the 
sub-arctic gyre. The gyres themselves are positioned as 
found in nature and are a result of the surface wind stress 
developed in the atmospheric model. The conditions in the 
northeast corner of the domain are similar to those found 
in the northeast Pacific. The warmer water is due to the 
convergence of the vertical shear current into the 



191 



northeastern boundary resulting in dovnwelling in that re- 
gion. The region of warm water and weak currents in the 
south central region of the domain is analogous to the water 
found in the Sargasso Sea area of the Atlantic. 

The zonal winds in the atmosphere develop a horizontal 
and vertical structure similar to that found in the real at- 
mosphere. A three cell meridional circulation was formed 
and slowly meandered several degrees of latitude north and 
south of its mean position. 

The second phase of the experiment consisted of a six 
year integration of the ocean-only model and a state of 
equilibrium was quickly reached. The upper level of the 
ocean had reached a state of almost steady temperature but 
the lower level still showed a very slow increase in the 
temperature. This was a good example of the extremely long 
time necessary for the vertical diffusion mechanism to bring 
the deep water into thermal equilibrium. 

Phase III of the experiment consisted of the integration 
of the atmosphere-only model equations. Energy balance 
studies of the atmosphere from this phase were made and com- 
pared with the same studies from phase I. The energy 
studies compared favorably with each other and with the 
existing knowledge of atmospheric energy exchanges. 

Data produced by the combined model and recorded at 
selected points vere spectrally analyzed yielding peaks in 
the spectra of both oceanic and atmospheric variables at 
periods on the order of one and two months. Cross-spectrum 
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analyses and time series analyses of several atmospheric and 
oceanic variables led to the following conclusion about the 
air-sea interactions. The longer period oscillations were 
present in the data and were attributed to air-sea interac- 
tion. An increase in the mean atmospheric flow preceded a 
rise in the barotropic flow in the ocean by about 3.5 days. 
This increased oceanic flow caused an increase in the advec- 
tion of heat and oscillations in the sea surface temperature 
about 7.6 days later. This change in sea surface tempera- 
ture led an increase in the vertical heat flux between atmos- 
phere and ocean and an adjustment of zonal mean atmospheric 
surface temperature. It is interesting to note that even 
with the limited model developed here, there are longer per- 
iod interactions between the atmosphere and ocean. If the 
model was modified so that more freedom were allowed, then 
interactions at different periods might occur. There are 
several ways in which this model might be changed for future 
studies. Some suggestions are listed below. 

The incorporation of continental land masses and the 
corresponding allowances for a zonal standing eddy into the 
model is a necessity. This effect would modify the heating 
in the model in that a more realistic longitudinal pattern 
of vertical heat flux would emerge. The effect of this 
change on the frequency of response in the model would be 
interesting to observe. The incorporation of a second har- 
monic into the atmospheric Fourier formulation would simu- 
late the continental effect. This would also cause a shift 
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in the position of the short period oscillation or add a 
new period of interaction. A study of an atmospheric model 
with a two-wave structure was carried out by Trumbov/er 
(1972), showing the feasibility of including a second v;ave 
in the atmospheric formulation. Incorporation of the second 
wave in the atmosphere would require an increase in comput- 
ing time but compared to the total time required for the 
experiment, it would be negligible. 

A change in the sensible heat term should be incorpora- 
ted in the model. In the term (T - T,), the quantity 

was held constant. This quantity should be divided into 
two parts, one that is constant and one that is a function 
of the surface wind speed. This would cause a portion of 
the vertical heat flux to be dependent on the wind speed 
and allow more freedom in the model resulting in other 
forms of interaction. 

An analysis of the energy budget of the ocean should be 
incorporated in any future calculations. The transforma- 
tion of energy between the vertical shear kinetic energy 
and the vertical mean kinetic energy and the dissipation of 
energy is an area well worth studying. 

The phase HI experiment should be carried out again 
with a sea surface temperature with an x dependence inclu- 
ded. A zonal variation in the heating would be present 
and the frequency at which the variables oscillated may be 
changed . 
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For atmospheric prediction for periods shorter than a 
few weeks, an interacting model may not be necessary. An 
atmosphere-only model with an accurat ely--observad air-sea 
temperature difference may be sufficient. However, longer 
range predictions on the monthly or seasonal time scale 
should include some kind of simultaneous ocean calculation. 

The atmosphere-ocean model was synchronously integrated 
for a 79 year simulation of the circulations of the two 
media. The interaction observed at periods of one and two 
months provides evidence of the importance of the effect of 
the ocean circulation on the transient behavior of the at- 
mospheric circulation, and vice versa. Other periods of 
interaction surely exist but the investigation of all of 
them is beyond the scope of this work. It is hoped that 
the results presented in this paper are sufficient to moti- 
vate further research in the important realm of large-scale 
air-sea interaction. 
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APPENDIX A 



Derivations of Equations (2.22) to (2.27) 



Equation (2.14) may be rewritten 



+ Ik X • V(v''-'j;j^j) + Ik X • 

3t 3t 

''(’"I't) + So ^ . 



A p 3x 9y 4 



where 



A(y>t) cos kx + B(y,t) sin kx, 
= F(y,t) + C(y,t) cos kx + D(y,t) sin kx , 



and 



■ '’S 1^4 I “« 
- PS Kl ''4 



(Al) 

(A2) 

(A3) 

(A4) 

(A5) 



An examination of quantities in Equations (A4) and (A5) 
shows 

^4 ~ 3y ’ 

^4 3x ’ 

and |VJ = [(||)2 + (|^)2] . (A6) 

Linear interpolation was used to obtain 

" f ^3 ~ J ^1 » 



and the relationships 






’^1 + ^3 
2 

2 
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♦l ■ ♦m + ''’t 



and 



Then 



^3 ^ 



(A8) 



^4 = ’I'm "^’^t • 

Apply Equations (A2) , (A3) and (A9) to (A6) to obtain (v/hen 

2k and higher order terms are neglected) 



3E ^ . 9F .3A . 9C. , 

u, = - T— + 2 -tt- - (t~ - 2 TT-;) cos kx 

4 9y 9y 9y ay 

,9B ^ 9D. . , 

- ^ 3?' 



(AlO) 



= k(B - 2D) cos kx - k(A - 2C) sin kx, 



and 

I'V. 



4 c 



M 4. irM _ 9M'>2 

3y " 8y'' 2^9y " 3y'’ ''3y 3y-^ 



+ (Ak - 2Ck)^ + (Bk - 2Dk)^] 

3E 3A , 3E 3C _ ^ 

^ 3y 3y “ ^3y 3y 3y 3y 

j „3C 3F. , , .„3E 3B , 3E 3D 

+ 8-^ COS kx + — 7 ^ 

dy dy dy dy dy dy 



, 9F 9B ^ o9D 9F. 

- 4-^ -TT — + 8-t : — -7^ — ) sxn kx 
dy dy dy dy 



1/2 



(All) 



(A12) 



Equation (A12) is simplified by making the assumption 
that contributions of the disturbances are much less than the 
contribution of the zonal wind to the magnitude of the sur- 
face wind. The cosine and sine terms in (A12) are therefore 
neglected and a new formulation for [w , | is obtained. 
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(A13) 



Iw, 



,dE Ir /3A 2 

+ 2 1 - 2-§7) 



+ (|^ - 2|^)^ + k^CA-2C)^ + k^(B-2D)^] 
9 y 9y 



1/2 



Equation (A4) is now 

I- I ir 3E , ^9F-, /9A „9C. , 

''x" ‘’'^D !''''« I 3y ^3y' ' ^3y ' 

- (|i - 2|£) sin kx 

9y 9y 

and Equation (A5) is 

T = pCjj K ^4 I {k(B-2D) cos kx - k(A-2C) sin kx 



9t 






9t^ 

9x 



The term [ 



9t 



9x 



9t 

*r 1, can be written 

9y 4 



9y U 



l_rhv I (2^ - -^)] 

3yU\V4l 3y^J 



, (A14) 



(A15) 



+ I^KVJ k^2C-A) - | 7 l|Wj( 2 |£ - |^)]j cos kx 
+ [I'Vj k^(2D-B) - fj-l IW 4 I (2|| - |2))j sin kx | 



(A16) 



When Equations (A2), (A3) and (A16) are substituted in 

Equation (Al) and terms of wave number 2k and higher are 
neglected, an equation of the form 

y^[G^(y,t) + G^ixy-y yt) cos kx + G^(x,y,t) sin kx] 

= Fj^(y,t) *f F^Cxjyjt) cos kx + F^Cx,y,t)sin kx (A17) 

is obtained. Coefficients of the terms independent of x, 
coefficients of the cosine terms, and coefficients of the 
sine terms are equated to yield the following three 
equations : 
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D 






3y‘ 



k^A] = kl 



^)E a^D a^E 



ay 



ay 



+ 

2 3y 



B 



9y^ ay 



i!l D - B - |£ D)] - s Bk 

ByS 9y 3y o 



A p 



pc 



•I 



IwJ k^A-2C)+ (2|^ 



’S”) 



a ^a^B 



tv' 



ay 



, - " - >‘t-3y 



3E a^A ap a^c , a^E 



8y- 



3y^ 



3y^ 



a 



ay 



C + k^(|^ A + C)] + e Ak 

3 3y 3y o 






3D 3B, 



pc„ J |\V^ I k^ (B-2D) 3- |^[ |w J (2|^ - 1^) : 



and 



a r3^E, k a y.a^B „3 ^a 

•^[ ] = T -^(A - B 

ay2 2 3y2 8y2 



+ c 



a^D ^a^c 



ay 



- D^) 

3y" 



+ 



A 



Ap os f? - If)]- 



Rewrite Equation (2.15) to obtain 



— (V2-y2)4)^ +|k X Vi>^ • 



dlpr 



+ lk X V 4,J . V(V^V -k 3 ^ ^ 



„ , 3t 3t 

- — £~ ^2 r — ^ - — —] - r- 

A p r^ 3 x 3 y -'2 '■ 



3t 3t 

J1 „ 






3x 3y'^2 ‘•3x ay"^4 

and examine the term y^KQ/2f . Substitute for Q to 



I 



u^kQ _ y^Kg 



2f 



2£^ - 2rA-r S'Sca-'-ir S ^ -(ff)<''i’=)]- 



, (A18) 



I . (A19) 



(A20) 



(A21) 

obtain 

(A22) 
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Define a constant, such that 



Cl = t“0.3 - c|^)A 



and the term is then 






2f 



2f Ap ^h'-'^sea 



R 



■(ip^ + C, ) ] . 



Substitute Equation (A3) for 4^^ and 

2f . 

[(F+C^)+C cos kx + D sin kx] 



y^KQ _ p^Kg 



H, 



2f 2fAp lisea R 
o o f 



is obtained . 

9t 9t 

V X 

The term -g/ A p 2 [ - g- --] ^ is then examined. 



that 



Tj " (P>;)2 02<"'3 



that can be broken up into its vector components 



■Xj " <“3 - “l> 



and 



T - (PK)^ (^^>2 (V3 - V^) 



The quantity roay be written as 

n 

( x ^)2 = ( PK >2 (^>2 ^ 



that leads to 

9t ^ 

= 2(pk)„ (^), ^ 

9y 2 p 2 g ^2 

The quantity may be written 






94^. 

9 X 



(A23) 



(A24) 



(A25) 

Recall 

(A26) 



(A27) 



(A28) 



(A29) 



(A30) 



(A31) 
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which leads to 





9Xy 

9x 


-2(pO, C^), 


(A32) 


and 


9t 

^ 9x^ 


- - - 2(P^)2(“>2 


(A33) 


Finally 


_S_ 
A p 


9t 9t 

2f ^ ^1 = 4 (dk) (^) V^il; 

^9x 9y ■'2 A p 2*^ p^2 


(A34) 



and can be v/ritten in its Fourier components as 

“7^ ^ (pK) ^ ^ “ (Ck^ “ ^)cos kx-(Dk^- ^~^)sin kx]. 

^ P 2 p 2 3^2 3y2 

(A35) 

Repeat the procedure previously used on Equation (Al) 
on Equation (A21) to obtain three equations. These are 






tl'TT ■ ’‘'sy 



3E 9^D . 9F 9^B 



9y 



I 7 B(k'- y^) - 



9 

9y' 



9y 
D - 



2 9y 



9E 



2 - D(k + y2) 



9 



9y 



9y 

B] - 3 Dk 



2 2f 

y KR 

2f A p 

O 



Hi(- -^C) + 4 (pk) 2(^)2(^ - Ck^) 

9y 



aV 



P 

9 r I.T7 



9C 



K|k (2C-A)- |-[lvvj (2gy 



9A 



9^D 






9y- 
9E 9^C 
' 3y^ 



9F 9^A . 9’f, . 

^ + A + 

9y <^„2 ^..3 



_ IL^ 

o 



9y 

2f 



9y 



9 ^F 

9y' 



4 ’ ^ 9y ay 
A(k" - y^) 

C] + 3 Ck 



)] 



] 



hyr "y- -r”> 



9^D 



j- S 

A p 



PS Ij'Sl f7'KSl(2|f - If)]] 



(A36) 



(A37) 
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and 



at "a 

9y 



- y"]F = j 



a , a^D 

ay^ 



•A - B 



a^c 

3y' 



+ - D^— ^] + |J^ T- fBC - AD] 

,2 -^,.2 oy 



ay' 



ay' 



- — H [T 

2f A p 1 sea 



2f 



■^(F d- C^)] 



+ _£_ 4(pk) (££) llF + -^ PC ^\\V I -2-^)1 
Ap p''2 . 2 A p P^D3y ^ I ' ''ay ay'* ^ ' 

dy 



This completes the derivation of the equations for 
dieting new values of the Fourier coefficients. 



(A38) 



pr e- 
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APPENDIX B 



Derivation of 0) and Equations, 



Applied at level two, Equation (2.12) is written as 

K 



*3 ■ ♦l 



^ ( 

9t A> p 



)+ Ik X • V ( 



^ ) + — t>)2= - f Q2* 

^ o o 



Substitute from Equation (2.13) which yields 

1^ + Ik X ♦ Vip^ - ff-E- ^2 = ^ Q 2 . 

o o 

Solve for which gives 



(B2) 



2f 



“2 0 A p 






8t 



3y 3x 



3x 3y 



it ^2 



(B3) 



Substitute Equation (2.16) and (2.17) into (B3) and 
neglect terms of wave number 2k and higher which gives 
2f 



“2 = 



w ) 1 rwT,8C ^3A, , , ,^3B . 3D, ^ 

FT7 - "^37^ 



or 



, 3F *^^2 , ,3C , ,3E^ 3F^. , , 

3T ■ 2 r •*■^37 ■ 

.,3d . , ,3E_ 3F. . , . , 



^2 = w + (u) ) cos kx + (w^) sin kx . 



(B4) 



The quantity (O 2 is composed of a mean component and a 
disturbance component. The mean component may be written 
f 



to = 



o ,, ,„3C ,.3A, . , ,^3B . 3D. 

oTJ - V' ‘‘^^'57 ■ V’ 



+ 



3F 

3t 



KQ; 

IT 



(B5) 
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and the disturbance part written as two components 

r 3C , ,3E_ 3F„. , 

w_ = — ^ k(-^D - -^B)] 

C a A p 9t 9y dy 



(B6) 



and 



“d “ 



11 ^ |i!> + k#c - |£a)] 

a A p '• 3t ^3y 3y 



(B7) 



The mean meridional motion component can be derived 
from the zonally averaged continuity equation. This is 
written 



r; X 

V • W 



3p 



= 0 . 



(B8) 



Applied at level one this yields 





(B9) 



and this is solved numerically to obtain v^ at any point in 
the model, after obtaining ^2 from (B5). 
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APPENDIX C 



Derivation of Energy and Energy Transformation Equations, 

To derive the energy relationships, the A, B, C, D, E 
and F equations [(2.22) to (2*27)] were multiplied by A 
through F, respectively. After further mathematical mani- 
pulation, the energy equations were obtained. 

Multiply the A equation by A to obtain 



and the B equation by B, and add to Equation (Cl) to get 




1 , 9^E^ 

] ~ IcA [ B 





- I^D) 3 - B ABk 
dy dy o 




(Cl) 







(C2) 
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Treat the C and D equations similarly and add to obtain 






- Ck - Cy^] + D^[ 



9 r3 h 



3 y^ 



- Dk^ - Dp^] = 



«y ■^„2 ^„2 



9 y 



9y' 






9 F,„ 9 ^B ^ 9 ^A, 



9y' 



- 

9 y' 



■^— ^(BC - AD) - (k^- y^) l^(BC - AD) 

9 y* ^ 5 ' 



P^B ||wy kh(2C^-AC) + (2D2-ED))-c|^[|\V^|(2|2 - |^> ] 



- lf>>| f? MpK) 2 (<lf) 2 (C( 



9 ^C 



9y- 



-Ck^) 



+ D( 



9 ^D 



9 y 



- Dk")]. 



(C3) 



Combine Equations (C2) and (C3) and integrate over the 



region to obtain 

y 



n 



a 9 f 9 a 

^ 9 Tt -7 

9y 



,9 B 



-Ak ]+B-^[- 

O t rs 2 

9y 



-Bk^]+c|^[^ 

Otr .2 

9y 



-Ck^-Cy^ ]+ 



p|ti— 

d t ^ 2 

9y 



Di^) 

9 y" 

+ 91 (a 11 d 

^ 9 y ^^,..2 



-Dk -Dy ] > dy= 



/ |“s"g 



9 ^B „ 9 ^A 



-B- 



+ C 



9 ^D 



9 y 



9y- 



- D 



9 ^A „ 9 ^C 



-B 



. 9 '"B, . ^ 29 F 



9y 



9 y 



tc )ty^-^(BC - AD) J 

9 y= 3 y^ 



+ f? pCB[(A-C)(hvJkhA-2C)+ |-l|,yj(2|| - M) 

+ (S-n)(l\vykhB- 2 D)+ |^l|ffj( 2|2 - 

2 2 1 

+ ^ 4(pk) (^) [C(-^ -Ck^)+D(^ -Dk^)] S dy. (C4) 

^ ^ 2 P 2 9 y^ 9 y 2 j 
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Perforin the multiplications and apply the boundary con- 



ditions to acquire 



71 ^'- 



|F(^^ _ j^ILa 

3y^ 3y^ 9y^ 3y^ 3y^ 3y^ 



- + C^)-y2|I(BC-AD) ]- 

3 y 2 3 y 2 3y Ap 



pCj^ j^(A-C)^|W^ lk^(A-2C) 



+ ' l7>' +(B-D)(|\Vjk^B-2D)+ |-(|wj(2|i 



- If'j)] - -m-y)l|dy 

(C5) 

Apply similar techniques to the E and F equations to 
obtain 



It/ u|f)^(|f) W hv = 



y 

71 



k(9E(Aiil _ gllA ^ (,3lB 



2 3y 3 y 2 



3y 3y' 



_ D— )+ —(A— -B-— -B— +C— )+y^|^(BC-AD) ] 
9y 'y-2 'i,.2 C5..2 :^„2 dy 



3y^ 9y^ 9y^ 3y' 

2 2f 



+ "ll’-sea- 4 PSl(>=-">V l'7 

^ o 

- If” - sf '■(P”2(“>2“-0’ I 



(C6) 
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Define 



- zonal potential energy, 

= eddy potential energy, 

= zonal kinetic energy, 

= eddy kinetic energy, 

C = conversion of P to P , 
p z e 

C = conversion of P to K , 
e e e 

C, = conversion of K to K , 
k e z ’ 

C = conversion of K to P , 
z z z 

= generation of P^ by diabatic processes, 
= dissipation of by internal friction, 

D^^ = dissipation of by surface friction, 

D^^ = dissipation of by internal friction 

and 

D = dissipation of by surface friction, 
z s z 



From Phillips (1956) and the basic definitions for ip. 



and used in this work (Equation 2.13) equations for the 



potential and kinetic energy components can be written as 

dy = dy, (C7) 



and 



/’ 



P 3 dy = f 



“V - I / 



CC + D ) dy. 



(C8) 

(C9) 



JK ^ 2 J 



•dy' 

2 2 2 2 2 
+k CA "+B +C +D )]dy 



(CIO) 



208 



The right hand sides of Equations (C5) and (C6) yield ex- 
pressions for energy generation, conversion and dissipation 
as follows: 



-/I 

‘■-7 



2Apf 



H, [T 
1 sea 



2f 
c 

R 



(F + C^) ]F \ dy, 



[| y^(BC- AD) ly] dy, 

9^A 



ceil) 

(C12) 



(k,9E ,.9^B 
9F 



+ C 



9^D 



+ (A 

9y 



9y^ 

9^D 



- D 



9y' 



^ei= - 



^s= - 




9y" 

9"A 

9y" 

. 2 , 



9y' 



- B 



9^C 



- D^^) 
9y" 

,9^B 



9y 



+ C-^)] Cdy, (C13) 

9y" ) 




D 



2 1 



and 



D = - 

2 S 




f- 4(p«)„(^)„[C(-^ -Ck^)+D(-^ -Dk^) ] (dy, (C14) 

2 p 2 I 

AF [(A-C) (|j?Jlc^A-2C)+ |y]|>vy(20 - M)]) 

+ (B-C) (^|\?^ I (B-2D)+ |j[ Kv^ I (2-|£ - |£) ])J dy , (C15) 

4(P>02(^)2 



(C16) 




£_ 

Ap 



pc, 






dy. (C17) 



The quasi-geos trophic first law of thermodynamics, 
Equation (2.8), can be manipulated in a similar manner to 
derive expressions for and . These are 

(C18) 



2 

and 
C 



■-/ 



2f 



Ap 



— CO F dy , 






D 



CC19) 
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APPENDIX D 



The difference form of the equation for the prediction 
of the vertical shear current, u\ analogous to Equation 
(3.8) is 



u: - u ' . 

ijk ij k 

2At 






_ Z i ^ fA \ 

^ijk " H ^ijk 



K = 1 



+ AV.f. + 

k j ijk ijk 



(Dl) 



where 



0^ 



ijk (Az). 



-dn). . 



AA 



[w . .+w . , ^ .+w . , , -+w . , ^ . ^), 1 



1 ijk 4 ij i+lj ij + 1 i+l,j-l k- 



(u),^ i. ~ •+'^•0.1 dw. .J^T+W^. 

'^2 



k--^ ^"ij ■ "i+lj ■ ”ij+l' "i+1, j+l^k+-^ ^'‘^k+-^^ij 



)..^i 



^ Ay(A.)^ + 



I 



-T-T — T — \ ( A z ) , A [(u. - u..)-(u.. - u. ) 

(Az) ^k i+lj ij i-lj 



+ (u.^ ^ - u..)-(u.. - u.. ,)], 

ij+1 ij ij ij-l'-'k 






-] 



(D2) 



and 



AVj^ = Ax Ay(Az)j^ 
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and 



AA = Ax Ay 



and 






- rl 






Ay(Az). 

8 






+ (u.. + U.. t) + (u.,T.,T + u ) 

IJ lJ-1 1+lJ+l 1+lJ 



* ‘“ij+1 “ij' k' 






Ay (Az) 



(u . . + u . . ^ ) 

ij iJ-1 



+ (u . ..+ U. .. -)+(u..,T + u.,) + (u. 

1-lJ lJ+1 13 ' 1-lJ+l 



+ “i-uMk' 



+ [-:7(u,.+ u..,^) 

'■2 xj iJ+1 



Ax (Az) 
8 



(v..,T + V. 1 .,,' 

lJ+1 1-lJ+l 



+ (v..+V. ,.)+(v.,T.+V,.)+(v.,T. T+V,. t) ),] 

ij i-lj i+lj ij i+lj-1 iJ-1 



+ “ij) 



Ax (Az) . 
8 



(v. . + V. T .) 

ij i-lj 



xj-1 1-13-1^ 1+lj ij' 1+I3-I ij-1 /k 

(D 3 ) 

The difference form of the equation of the meridional 
component of the vertical shear current, v*, analogous to 
Equation ( 3 . 9 ) is 



j-jk ijk 

2At 



= 



ijk 



1 ^ £ 

- I E ’‘ijk 



k=l 
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■AV f . [a(u') + (1+a) 

k j ' ^ ' ijk ^ ' X 2 k 



(D 4 ) 



where 

ijk 



AA, 



(Az). 



-do),- 41, - —[ (W 4 4+"44.1 4+W.: 4J.1+W -.-4.1 44-1 )u_i. 



^'ijk 4 " i j ' " i+1 j ’ " i j +1 ' " i+1 j +1 '' k-i 



( v) , 1 - (w . .+w ... . +w . . , - +w (v),,!].. 

k-j ij i+lj iJ+1 i+lj+l'k+y ^k+y-’ij 



- ^ Ax(Az), t[ (P4 44.1 +P,- 4.1 44.1 ) - (P4 4+P44.1 4) i 



k 2 ' ij+1 i+lj+1 



ij i+lj k 



i 



+ 7-7 — r — ^ (Az) , A [(v. ^.-v..) - (v.. - V. ^.) 

(Az), I 'k K'-^ i+lj 12' ij i-lj 



+ (v. - V..) - (v.. - V.. ,)], 

ij+1 IJ ^ ^ IJ iJ-1 k 



Idv - V ) k(v - V ) i£-l 

+ AA[ 77 -ir 7-^ , ■ ■ ] 



(Az)k_l, 



and 



(I^) . ;zn~(v. .+v .). ' 

2 ijk j 2 xj x+lj k 






Ay ( Az ) 
8 



(D 5 ) 



(u..^. + u.,T. t) 

i+lj i+lj-1 



+ (u..+ u.. ^) + (u.,^.,-) + (u..^+ u..),] 

IJ xj-1 1 +lj+l xj+1 XJ Jk 






Ay ( Az ) 
8 



(u,.+ U.. t)+(u. t.+U. t. -) 

IJ XJ -1 x-lj 1-lj-l 



+ (u.. ^+ U.,) + (u. T.,,+ u. ,.)),] 

XJ +1 XJ 1-lj+l 1-lj Ik 






Ax (Az ) 
8 



(''ij+i'" '’i-ij+i) 



+ (v..+ v. -,.) + (v. .,.+ V..)+(v. -. ^+V.. i)|ll 

XJ x-lj-^ 1+lj IJ x+lj-1 xj-lMk-* 



+ ■"''ij)' 



Ax (Az) . 
8 



(v , . + V . T . ) 
XJ X-ll^* 



Ij 



+(v,. 1+v, ,, n)+(v.-,+v..)+(v.,T. , +v .. ))]>.(D6 ) 

xj-1 i-lj-1^ x+lj ij' i+lj-1 ij-l^yk-^ 
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At each time level H y Equations (Dl) and (BA) are solved 

jl+1 

together for the unknown shear currents, . and 

, £+1 
^ijk • 

The difference form of the vorticity equation analogous 
to Equation (3.12) is 

i +<^3)1-1, J-l) - 

+ — ((!) +(I) ^^-3 ^ 



^ -u . . ) 



(D7) 



where 


km 


(I ) ^ 

^ 3^ij 


E 

k=l 



- ("ij- + 



(it-1) 



(D8) 



<^4>iJ - I E \l><''i +13 - ''ij) 

k=l 



- (v.,-v. , . ) + (v . . ^ -V . . ) - (v . . -V . . ,)]5^ 

IJ x-lj iJ + 1 IJ IJ ij-1 k 



+ AA[ 



^^ ^k-l~ ^khi _ ^^^k~ ^k+l^ii j^-1 



(AOk_i 

2 






and 



(D9) 



IJ 2 ^l+lj+l ^x-lj+l ^x+lj-l ^x-lj-1 ^XJ 



(DIO) 



Equation (D7) is solved by relaxation techniques to 
5/+1 



obtain Ui . . 

ij 
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The difference form of the thermodynamic equation analo- 



gous to Equation (3.16) is 



_ (ji+i) „ a-i). 

^ i.ik ~ "i.jk 

2At 






+ w. + w. 1 (T), 1 - (w..+V7 ^.+W.^ ^ 

ij+1 1+lj + l k— 2 i+lj iJ+1 

+ w. (T), _^1]^ + (Az), A„[ (T. , - . - T. .) 

i+ij+1 k+Y ^ ^ ^+ij 



- (T.. - T. ,.) + (T...,- T..)-(T.. - T.. J], 

ij i-lj iJ+1 ij ij iJ-1 k 



(ii- 1 ) 



+ AA[-7- v ^ 7T^^ V • -■] • • ^ ^ 

(Az)j^_l_ 



k+T 



(Az) 



k 



At 



- 1 



+ 

where 



(I,).., =-!f[(T.,+ T..,.), • (u..+ u.. -), 

5 ijk ij i+lj k ij ij-1 k 



- (T . , .+ T . . ), 

i-lj ij k 



(u. ^,+ u. ^ . t), 

i-lj i-lj-l^k 



+ (T, .+ T . . . , 

ij xj + l'k 



(v. .+ V. , .), 
ij i-lj k 



an d 
(Az). 



At '“‘^C 
vious ly . 






(6T) is the convective adjustment discussed pre 



(BID 
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